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CAPÍTULO 1 


FUNDAMENTOS 
EM MODELOS 
DINÂMICOS 


1.1 SISTEMA 


O conceito de sistema pode ser definido de diferentes formas. Em 
controle de processos, denota-se como um objeto ou uma coleção de 
objetos que realiza certo objetivo e cujas propriedades pretende-se es- 
tudar. Alguns exemplos são: sistema de fabricação de papel ou cerá- 
mica, planta solar, circuito elétrico, servomecanismo de posição, sistema 
biológico ou econômico, manipulador robótico, reator, coluna de desti- 
lação, laminador, trocador de calor, refinaria, entre outros (AGUIRRE, 
2007; OGATA, 2010). 


A figura (1.1) ilustra os principais elementos de um sistema de 
controle. 


EN 
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Figura 1.1 = Componentes de um sistema de controle 
meto | CH | E 
5 
entrada (causa) H() saída (efeito) 


sistema 


Os problemas de controle associados à estrutura da figura (1,1) são: 

* Análise: é conhecida a entrada, u(.), o sistema, h(.), e deve-se 
obter a saída, y(.); 

* Projeto: é conhecido o sistema, h(.), a saída desejada, y() e 
deve-se obter a entrada, u(.), para proporcionar tal saída; 

* Identificação: é conhecida a entrada, u(.), a saída, y(.), e deve 
se obter o sistema, h(.), de modo que j(.) —» y(.), onde y()éa 
saída do sistema real (medida) e y(.) é a saída estimada. 


1.2 MODELAGEM E IDENTIFICAÇÃO 


Entende-se por modelagem e identificação a determinação do 
modelo matemático de um sistema real em estudo representando 
suas características de forma adequada para uma utilização particular 
(detecção de falhas, otimização, controle, entre outros). Os procedi- 
mentos envolvidos na elaboração de modelos matemáticos estão mos- 


trados na figura (1.2) e são usualmente classificados em (LJUNG; GLAD, | 


1994; BOBÁL et al., 2005): 
* Análise físico-matemática: baseia-se nas leis da física que 

caracterizam um sistema particular, tais como as leis de com 

servação de massa, energia e momento (modelagem fenomé 

nológica); 

Análise experimental: 


baseia-se nas medidas ou observações 
do sistema (modelage: 


mM experimental), 


EN 
K 


b 
& 
h 
t 


SEE à 


fz 
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Figura 1.2 = Princípios para a construção de um modelo matemático 


Estes procedimentos propiciam a obtenção de modelos que 


aço a dinâmica do sistema (processo ou planta). Para fins de 
- trole de processos, não se pretende encontrar um modelo ma- 
as temático exato, mas um modelo adequado para uma determinada 
eo pro Na prática, utiliza-se a hipótese básica de que, para a elabo- 
mo” ão “se processos reais, em geral, não necessitam obrigato- 
riamente de modelos complexos (HANG; CHIN, 1991; GESSING, 1996; 

LJUNG, 1999), 
quaté “O modelo de um sistema é uma equação matemática utilizada para 


tesponder a questões sobre o sistema sem a realização de experimen- 
' de um modelo pode-se calcular ou decidir como o sistema 
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16 Ra. 
ições operacionais). À utili 
-se sob determinadas cond é à 
modelo aims simulação do sistema quo ipi 
baixo custo e seguro para experimentar ed a pede Er , a 
dade (adequação) dos resultados de simulaç alidad 
do modelo matemático do sistema. 
Alguns dos diferentes propósitos para 
matemáticos em automação industrial são: 
* Predição: é uma tentativa de predizer os estados futuros dy 
sistema (comportamento dinâmico) e está limitada à 
do modelo e aos efeitos das perturbações atuantes (prese 


no sistema; 


= Análise e projeto de sistemas de controle: proporciona um 
vasto campo para aplicação em modelagem e identificação 
na sintonia de controladores clássicos, síntese de algoritmos 
de controle adaptativos e preditivos, dimensionamento de 
equipamentos e estimação do estado de variáveis não mensu- 
ráveis; por exemplo, a estimação da velocidade a partir da po- 
sição é uma medida indireta; 


* Supervisão: emprega a simulação, com base no modelo mate 
mático, para a avaliação das características operacionais do 
sistema, para o projeto de engenharia e para o treinamento de 
operadores. Muitas vezes é também utilizado na detecção 


erros, no diagnóstico de malhas de controle e no au de 
operadores em atividades v 


programa de simulação p 


comport 


a utilização de Modelos 


e a a 
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1.3 DESCRIÇÃO DE SISTEMAS: MODELOS 
CONTÍNUOS 


A representação de sistemas pode ser realizada pela função de 
transferência, resposta impulsiva ou através das equações de estados. À 
seguir, apresentam-se de forma resumida as três abordagens utilizadas 
para descrever sistemas dinâmicos contínuos (DORE; BISHOP, 1995; 
e JACQUOT, 1995; SEBORG et al., 2004; OGATA, 2010). 


a) Função de transferência contínua 

A função de transferência H(s) de um sistema é a relação entre a 
transformada de Laplace da saída, Y(s), pela transformada de Laplace da 
entrada, U(s), ou seja, 


HM E 


é uma razão de dois polinômios em s e está representada por 


HP 4.2 
E 


onde 


B(s)= »by FO A(s)= Dias 
Fo i=0 (1.3) 


sendo n a ordem do sistema (polos), m o número de zeros, a,=1en 2m 
(sistema causal). Os elementos (b, a, n) com je [0,m] e ic [0,n] são des- 
conhecidos e devem ser determinados analiticamente ou estimados. 
Algumas vezes assume-se que o valor de n é conhecido a priori. 


is o b) Resposta impulsiva contínua 
e Em um sistema linear invariante no tempo, os sinais de entrada, 
Ps uí(t), e saída, y(t), estão relacionados pela integral de convolução dada por 


ame = 
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y(0= [ b(t- DutDdr 


onde h(t) é a resposta impulsiva do sistema. 


A resposta impulsiva está relacionada com a função de t 


rência por E ] 
(o =L'[H(S) (14 


onde L“ é a transformada inversa de Laplace e t é o tempo contínuo, 


c) Equações de estados contínuas 
A função de transferência está relacionada com a representação de 
estados, no caso monovariável, de acordo com a seguinte relação: 


H(s)=c(sI-A)?b (3 


onde À é a matriz do sistema (nxn), b é o vetor de entrada (nx1), c é o vetor 
de saída (1xn) e Té a matriz identidade (nam). As equações de estados são 


x(t)= Ax(t)+ bu(t) 
y(D=cx(t) 


sendo x(t) o vetor de estados (nx1) 
do sistema. é 


(1.6) 


u(t) o sinal de entrada e y(t) a saída 
1.3.1 EXEMPLOS DE M 
DA FÍSICA ODELAGEM POR LEIS 


Representaçã 
ção matemáti i 
Equsisies mática dos Sistemas elétrico 


Considere a fi 
q. SR - gura (1.8) ue j 
Srcuito elétrico RLC e um Pis Pomba o Rá imples de um 
Sa-mo. amortecedor. 
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Figura 1.3 = Sistemas elétrico e mecânico 


u(t c y(t) 


K 


Sistema massa-mola-amortecedor 


Com base na lei das tensões de Kirchhoff e condições iniciais nulas, 
tem-se a equação que descreve a dinâmica do sistema elétrico, ou seja, 
u(=Rico+ LEOA +Yoficoar ; o=CVO am 


onde R é a resistência, L é a indutância, C é a capacitância, u(t) é a ten- 
são de entrada e y(t) é a tensão de saída. 


“Por outro lado, aplicando a segunda lei de Newton no sistema 
r rtecedor, obtém-se 


10 = e) + BOA mA, (1.8) 
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ij de atrito, M é à 
onde K é a constante da mola, B éo ea , Masga, 
ft) é a força (entrada) e x() é 0 deslocame 


destes sistemas (modelag, é 
EraGSÇT de Laplace são dadas Por 


As representações matemátii 
nomenológica) em termos das trans 
1 x() 1 


+, F(s) Ms? +Bs+K 


U(s) LCs? +RCs+1 


As funções de transferência obtidas com as Pair ciates ss físicas 
em termos das constantes R, L, Cou M, K, B são casos especiais dos 
modelos entrada-saída. Em geral, os modelos avaliados, nestes casos, são 
denominados modelos paramétricos e representam uma dada estrutura 
onde os parâmetros são algumas vezes desconhecidos e devem ser 
estimados. Por outro lado, para avaliar os parâmetros (valores numéricos) 
através da modelagem por leis da física, necessita-se conhecer as con- 
dições internas e externas bem como a estrutura física do sistema. 


(1,9) 


Exemplo 1.1 - Seja o sistema massa-mola-amortecedor da figura (1.3), 
Admitir que os parâmetros físicos do sistema são: M = 2 kg, B=4 Ns/m, 
K = 12 N/m. A resposta temporal da saída (deslocamento), para uma 
entrada F(s) = 1 N, está ilustrada na figura (1.4) com comportamento 


subamortecido e a Programação com o software científico Matlab, da 
empresa MathWorks, na tabela (11). 


Conforme mostra a figura (1.5), se 
mentado para B = 10 Ns/m, 
oscilatória e lenta (compo; 


1 9 coeficiente de atrito é incre- 
então a dinâmica da saída torna-se não 
rtamento sobreamortecido). 


amortecedor 
“ir parâmetros 


& Avali 
RR Fesposta 


er 
Bina polos de sigam 
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Figura 1,4 = Caracteristica temporal do sistema massa-mola-amortecedor 
mf 012 
7 


01 
008 
E 


0.06 


0.04 


0.02 


o 1 2 


Figura 1,5 - Saída do sistema massa-mola-amortecedor para 8 = 10 Ns/m 


D.09 


0.08 


rerê 0.07 
o 0.06 
0.05 
0.04 
| 003 
2 25 3 
A 
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Representação matemática do sistema motor DC 


Seja o motor DC (Direct Current) controlado por armadura a 
forme ilustra a figura (1.6). 


DC controlado por armadura 


2 otor 
Figura 1.6 - Diagrama simplificado de um mot 


iF constante 


As equações diferenciais que determinam o comportamento diná- 
mico do motor DC são expressas na forma 


100 p= ; TaD=Ki) 
(1.10) 


e(D=Kwít) ; Lo +Ri(t) + e(=v(t) 


Os parâmetros e as variáveis de interesse são: R resistência de 
armadura), L (indutância de armadura), i(t) (corrente de pomee io 
ps de ii v(t) (tensão de armadura — entrada), e.(t) (força 
contraeletromotriz), w(t, , 
fornecido pelo Ee, to do oral TO (A 


à carga referida ao eixo do motor), e 
A relação W(s)/V(s) é dada por 


Wís) E K 
“O ++ RITR E “a 
ap 


t4, 


Elst 
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A representação matemática e o comportamento dinâmico do motor 
DC são conhecidos se os parâmetros são fornecidos pelo fabricante (es- 


pecificações de catálogo). 
Representação matemática do sistema monotanque 
Considere o sistema monotanque conforme a figura (1.7), onde 
A é a área do tanque (m”), a é a área do tubo de saída (m?), hi(t) é o nível 
do líquido no tanque (m), u(t) é a vazão de entrada (m/s) e O Ea vizão 
de saída (m/s). 


Figura 1.7 - Processo de nível 


q 
a; ——s 


Para avaliar suas características operacionais, deve-se obter uma 
equação diferencial não linear para o sistema (calcular um modelo ma- 
temático relacionando as variáveis de entrada — u(t) e saída — hr(t)). 

Assim, utiliza-se a lei de Bernoulli que descreve a relação entre a 

de vazão da saída (m/s) e o nível do líquido no tanque, isto é, 


E v(D=/2gh, (O) (1.12) 
aceleração da gravidade. 


a 
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a E tg, 


A equação que relaciona a vazão da saída, q(t), e a velocidade á 
e 


vazão da saída, v(t), é dada por 
q(t) av(t) ty 


O volume de líquido no tanque, em um instante t, é calculado " 


Volume = Ah, (t) (m”) (119) 


e modifica-se de acordo com a diferença entre o fluxo de entrada e Saída 
(denominado balanço de massa), isto é, 


dah (O , 
REA q(t) + u(t) (149 


Substituindo-se as equações (1.12) e (1.13) na equação (1.15), 
obtém-se uma equação diferencial não linear para avaliar o compor. 
tamento do nível do líquido, ou seja, 


dh(t) -ay2g 1 
>= Vh(D+—ut(t 
Er TR vhi(t) auto) (116) 

Pelo conhecimento da dimensão dos diferentes elementos que 
compõem o sistema monotanque, pode-se avaliar sua dinâmica pela 


equação (1.16); por exemplo, é possível determinar o nível h(t) quando 
o fluxo de entrada u(t) é conhecido. 


O fluxo de saída é calculado por 


a()=a/2g./h, (1) (um 


14 DESCRIÇÃO DE SISTEMAS: MODELOS 
DISCRETOS 


A seguir, apresentam-se de fo resumida abordagem” 

4 rma 
utilizadas para descrever 5 “= 
Matlab, pode-se elaborar e 
(comando c2d) e vice-vers 


Bim 
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das equações que descrevem as transformações retangular ou trapezoi- 
dal, determinar os correspondentes modelos contínuos e discretos. 


a) Função de transferência discreta 
A função de transferência discreta é a relação entre a transfor- 
mada-z da saída, Y(z), pela transformada-z da entrada, U(z), isto é, 


H(a)=" Ya) (18) 
A relação H(z) é uma razão de dois polinômios em z e é repre- 
sentada por 
d 
H(g)=B SA 
(z) na(2) (1.19) 
onde 
Bi(7)=>biz' ; A(z)=Daiz' (1.20) 
Fo 1=0 


sendon >m, aí =1 en a ordem do sistema. Os elementos /b/,a/,n/ 
com je [0,m] e ie [0,n] são desconhecidos e devem ser determinados 
através de uma modelagem matemática do sistema ou através da iden- 
tificação. Algumas vezes assume-se que o valor de n é conhecido a priori. 


b) Resposta impulsiva discreta 
A resposta impulsiva está relacionada com a função de transfe- 
rência por 
h()=2"[H(2)] (1.21) 
onde Z4 é a transformada-z inversa e t é o tempo discreto. 


As amostras da resposta impulsiva, h(t), estão relacionadas com as 
amostras da resposta ao degrau, s(t), pelas seguintes equações: 


no =LIst=s(t=D] ; s0=hO) 
T 150 


Busto 
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c) Equações de estados discretas 


i a com a represe; 
A função de transferência está relacionada cc p ntação de 


estados discreta, no caso monovariável, por 
dylud 
H(g)=e(Zl= A") b ti 


onde 4º é a matriz do sistema (nxn), b! é o vetor de entrada (nx1), cé, 
vetor de saída (1xn) e I é a matriz identidade (nxn). As equações de esta 
dos na forma discreta são 

x(t+)= Atx(0)+b'u(t) 


YO =ex() sa 


sendo x(t) o vetor de estados (nx1), u(t) a sequência de entrada e y(t) 
saída (medidas especificadas a cada período de amostragem). 


Na avaliação de modelos matemáticos discretos, deve-se selecionar 
um período de amostragem, T, (sampling period), para cada aplicação par- 
ticular, de acordo com uma das seguintes relações (JACQUOT, 1995; 
SEBORG et al, 2004; ASTRÓM; WITTENMARK, 2011): 


ta/T=5..15 ; T=7/10 (1.24) 


onde tssx é o tempo que à resposta do sistema leva para alcançar 95% 


do valor de regime permanente á 
do sistema. aii € 7 é a constante de tempo dominante 


1.4.1 EXEMPLO DE MODE 
LAGE 
Seja o circuito RC da fi e 


é a tensão de entrada. Bura (1,8) onde (t) é a tensão de saída eult) 
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Figura 1.8 = Circuito elétrico RC 


R 


u(t) cl yt(t 


A equação diferencial linear de primeira ordem que representa a 


dinâmica do sistema é dada por 
dy(t 
ReHO., y(t)=u(t) (25 
dt 
A partir da aproximação numérica da derivada (Fadali e Visioli, 
2009) dada pela equação (1.26) 
dy(t)  y(D=y(t=1) 

4 des D: Ea 
obtém-se a função de transferência discreta para o circuito elétrico RC 
de acordo com 

0 (1.27) 


U(z) (RC+T)z-RC 


Selecionando uma tensão de entrada para o circuito elétrico RC 

(que pode ser obtida por uma fonte variável DC) e conhecendo-se os 

numéricos de R, C e T,, as operações de carregamento e descar- 
Tegamento no capacitor podem ser avaliadas. 

Outro exemplo é o tanque para aquecimento de um fluido descrito 


na figura (1.9). 
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cesso térmico de malha aberta 


pa 


Figura 1.9 = Prot 


amplificador 


O tanque contém R ecimen! 
cai um dispositivo de i 
meg gonna saio NA a A 


(MC/ua) = E p | 
5.8 min e 9. = 1.5 min, 


co 
rrespondente 
desprezando-s de transferência 
Eid 
EP x 1/05 +1), é caraçto A temperatura ta de segunda 
“aracterizada pela te e peteão os 
E 7 | 
Co GR 
9 ASsrnas as 
8s+1) 
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Admitindo a presença do segurador de ordem zero (associado nó 
conversor digital-analógico) em série com o modelo contínuo do pro- 
cesso térmico, tem-se a representação discreta de segunda ordem (pe- 
ríodo de amostragem de 1 s), ou seja, 


T(z) 0.00832+0.0065 
Q(z) 2º-1.45212+0.4819 


que pode ser obtida pelos seguintes comandos do Matlab: 
gps= tf(0.496,conv([1.5 1],[15.8 1])) ; gp z=c2d(gp s,1,zoh) 


1.4.2 EXEMPLO DE MODELAGEM POR ANÁLISE 
EXPERIMENTAL 


Considere o processo discreto de primeira ordem caracterizado 
pela seguinte equação a diferenças: 


y(t+1) = 0'y(t) + u(t) (1.28) 


onde 9” é o parâmetro desconhecido e y(t) = O (W't < 0). Admita o se- 
guinte modelo para estimação: 


p(t+1) = Oy(t) + Ult) (1.29) 


sendo Ô a estimativa de 0" e (t+ 1) a predição da saída ou o valor da 
saída no instante (t + 1) baseado no parâmetro estimado ô. 
Seja a função custo dos mínimos quadrados dada por 


Jt)= EO) (1.30) 
20 


e(t) = y(t)— y(t) (1.31) 
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30 — — E 
(1.28) e (1.29) na equação (1.31), Obtémgg, 


Substituindo-se as equações edidas de entrada e saída (atual é 


m 
erro de estimação em função cats 
anteriores) do processo de aco! 


«0=0y(t-D-Oyt-D=y(0-utt=D-y(t=1) qa 


i é conhecido o parâmetro |" | A 
ra certo sistema que não : 
as Ni é entrada e saída no intervalo da perriapeçã su 
Assim, obtém-se o estimador dos mínimos quadrados p: encia 
ção de J(t) em relação a O, OJ(t)/00, resultando 
t=1 
Eytobye+1)-u(o)] 


Gg —— (133) 


yo 
k=0 


Tabela 1,2 - Programação em Matlab do estimador da equação [1.33] 


timador do parâmetro da planta da equação (1.28) 
somal=0; somaZ=0; 

u(1)=0; u(2:N)=1; 

for k=2:N % Obter medidas 
YIK)=0.9*y(k-1)+u(k-1); 

end 


1 $ Estimar Parâmetro 
Somalt+y(k)* (y(k+1)-u (x); 
SomaZ=somaZ+y (k) au 


Íeta (k)=somal/soma2; 


teta(k),+ 


ai e E 
a) k', linewidtn,2); 


l('tet, *label ('amostras') ; — em 


O tempo discreto éui 


9 instante +, À estimação do pa 
Pela equação (1,33) minimiza (1.30) pela suposição de 
PR Pri (1.28), De acordo ps 
9 e necessita de todas Pipes prnindo e saída 
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dá experimentação (procedimento off-line). A solução da equação (1.33) 
(estimativa de 9) pode ser obtida com um software numérico, tal como 
o Matlab. A tabela (1.2) ilustra a programação em Matlab da equação 
(1.33), onde se admite 0" = 0.9, uma sequência degrau unitário e o tem- 


po do experimento de 30 amostras. 


1.5 PROBLEMAS 


1) Considerando a constante de tempo mecânica (tm = J/B) maior que a 
constante de tempo elétrica (t. = L/R) de um motor DC controlado por 
armadura, obter a equação diferencial de primeira ordem em termos da 
velocidade e tensão de armadura. 


2) De acordo com as especificações de um fabricante que comercializa 
um motor DC, tem-se a seguinte parametrização: R = 20 (2, J = 2 Nms”/rad, 
B= 0.1 Nms/rad, K, = 1 Nm/A, Ky = 3 Vs/rad. Avaliar o comportamento do 
processo por simulação para as seguintes condições: 

a) Admitir condições iniciais nulas, 7, >> 7; e aplicar um degrau unitário 
na tensão de armadura para observar o comportamento dinâmico da 
velocidade; 

b) Observar as dinâmicas temporais analógica e digital (período de amos- 
tragem de 0.1 s) sob as mesmas condições de operação do item (a). 


3) Avaliar por simulação o comportamento do sistema de nível, equação 
(1.16), quando uít) = 1, t >0 (degrau unitário) e ho = O. Observar 
também a característica temporal para u(t) = O e hr = 0.25 m. Admitir 
A=le a/2g =1, Adicionalmente, considerar o ponto de operação (to, hi) 
£ mostrar que o modelo linearizado é de primeira ordem, relacionando 
Hds)/U(s). 


4) Simular os comportamentos de carregamento e descarregamento da 
no capacitor do circuito RC a partir dos modelos contínuo e 
- Admitir R=1MQ, C=1uPeT,;=01e1s. Se ão de 
s=0. . Seja a tensão 
entrada do circuito de acordo com 
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1 0<t<10 
u(t)= ô t<10 e t>10 


em Matlab da tabela (1.2) para estimação q, 


5) Implementar o programa linear discreto de primeira ordem ra. 


parâmetro 0 relativo ao processo 
ntado por 

Rouge y(t) = Oy(t= 1) + u(t—D) 

onde 9= 0.75, u(t) é uma sequência degrau unitário discreto e sendo 29 

amostras o tempo total da experimentação. 


6) A tabela mostra os valores medidos da resposta impulsiva de um sis. 


Para fins de aplicação de um controle digital direto, qual deve ser o pe 
ríodo de amostragem do sistema discreto de malha fechada? 


7) Considere o problema da estimação de um sinal constante a partir 
das medidas ruidosas, isto é, y(t)=s(t )+eft), sendo e(t) o ruído do sis- 
tema. À estimativa natural do sinal é calculada por 


a) Reescrever a estimativa de modo ) 
quando uma nova medida y(t) está di te ado eia RO 
e essência de um algoritmo de esti gi 
Implementar o resultado do i 
te; 
à convergência da estimativa, d nha o ae Mt a 
mando rand na geração de por * 100 amostras e utilizar oc” 
€) Comparar os resultados n : 
Passa-baixa digital de Primeira or “era implementação de um filtro 
Medida y( o 
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8) Obter a função de transferência H(z) = Y(z)/U(z) e a resposta impul- 
siva, h(t), associada à seguinte equação a diferenças: 


y(D=1.4y(t-1)- 0.4y(t-2)+u(t-1) 


9) Determinar e esboçar a resposta temporal do sistema discreto caracte- 
rizado por 

y(D=y(t-1)-0.09y(t-2) +u(t—1) —0.82u(t—2) 
Admitir como entrada uma sequência degrau unitário e condições ini- 


ciais nulas. Qual o polo dominante e o valor da saída em regime per- 
manente do sistema? É o sistema estável e de fase mínima? Por quê? 
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NOÇÕES 
BÁSICAS SOBRE 
IDENTIFICAÇÃO 


2.1 CONCEPÇÕES PARA IDENTIFICAÇÃO 


A identificação de sistemas é tratada, muitas vezes, como um 
problema de otimização que envolve algumas medidas para adequa- 
ção de modelos candidatos a representar um processo real. A seleção de 
modelos matemáticos e o ajuste dos parâmetros são influenciados por 
diversos fatores, entre os quais: (a) conhecimento a priori do sistema (li- 
nearidade, grau de não linearidade, atraso de transporte); (b) proprie- 
dades do modelo do sistema identificado (complexidade); (c) seleção da 
medida do erro a ser minimizado; (d) presença de ruídos (JOHANSSON, 
1993). A identificação de sistemas é um claro exercício que envolve 
múltiplos e conflitantes objetivos, tipicamente complexidade do mode- 
lo, critérios de desempenho e validação, que influenciam a seleção das 
estruturas do modelo matemático, A noção de um “bom modelo” é subje- 
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UPE tiva e erro seja uma regra relevante q 
m que a tenta a (ISERMANN, 1980; LJUNG, 199 


genharia em identificação de proces a ordem do modelo tão p.; 


sm para manter É - 
Várias razões existem ENO podem ser introduzidos ; 
ção com O princípio da Parcimônia 


tiva, fazendo co) 


op, Ci de fo pro 
combinar à necessita — rérios de informação Bayesiana, de Ali 
Estes critérios, tais como critérl lada valha residual e a E 
ou minimum description ver 1990). A meta do algoritmo E E, 
do modelo (HABER; UNBE maximização de um critério de dese, 
RC trt de acerto ou representação do siste 
penho que está associado à taxa de acert( EEE Ma 
avaliado. Se todas as restrições e condições forem aten s, o Modelo 
encontrado pode ser aceito. Caso contrário, se ia das condições im- 
postas é violada, todos os procedimentos de identificação, estimação de 
parâmetros e diagnóstico do modelo devem ser reavaliados até que um 
modelo apropriado seja encontrado (LJUNG, 1999). 

Muitos dos métodos de identificação, tais como os baseados nos 
estimadores dos mínimos quadrados ou máxima verossimilhança, são, 
em essência, técnicas de busca local guiada por gradiente e necessitam 
de um espaço de busca regular ou um índice de desempenho dife 


7: a em de algumas desvantagens, tais como: 
(a) alguma informação inicial dos Parâmetros do pisçr é necessária 


o 
LJUNG, 1996; BOBÁL et al, 2005 


Fundamentalme , 
minação de um modelo mater cação de Sistemas consiste na det” 
pah do sistema, caracterizado dios *epresente os aspectos 
e 7 i É 
pis aro eis Telacionados atravéo de ação doe einaia de entinl 
ego screta (ISERM, ELAÇ) * uma função de tran: 
tratamento a a industriais, o N, 1985; LJUNG, 1999) 
medidas (procedi nto elo pode ser obtido a parti ã 
estatístico, filtragem de 
» filtragem 
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coletadas a partir da realização de um experimento. O modelo matemá- 
tico final é uma forma do conhecimento da relação existente entre os 
sinais de entrada e saída, caracterizada no processo físico pela função de 
transferência. A figura (2.1) ilustra a composição básica em blocos de 
uma tarefa de identificação. 


Figura 2,1 - Procedimento para identificação de processos 


incertezas 
entrada | saída 
processo dinâmico —(4) 
técnica de 
ado identificação 


modelo matemático 


de ur A área de identificação tem tido considerável interesse nos 

|. últimos anos para fins de predição, supervisão, diagnóstico e controle. 

est Observa-se a aplicação em diversos campos da engenharia, tais como 

o sistemas químicos, térmicos, biomédicos, mecânicos, econômicos, elé- 

tricos, aeronáuticos e outros mais. Em uma tarefa de identificação, 

diferentes procedimentos para geração do sinal de entrada, medição da 
saída e armazenamento dos dados são utilizados, entre os quais: 


il * Identificação de um processo pelo teste da resposta ao degrau: 
39 9 processo é submetido a uma mudança na entrada do tipo de- 
grau e o armazenamento das medidas é realizado por algum 
registrador. Com a curva de reação do processo é possível 
aplicar diferentes técnicas gráficas, numéricas ou computa- 
cionais para modelar o sistema controlado por funções de 
transferência de primeira ou segunda ordem (processos com 
modelos de ordem reduzida), conforme ilustra a figura (2.2). 
SA O teste da resposta ao degrau só tem validade para processos 
a gt lineares ou para processos não lineares que sejam aproxima- 
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na vizinhança do ponto de PP tata te 
E or natureza, não permite a estima, 
já que o sinal degrau tem uma 


damente lineare: 
da resposta ao di 
de modelos de or 
pobre composição e! 


legrau, por D 
dem superior, 
m frequência. 


Madelagem a partir de uma mudança na entrada 
Modelagem a pº 


Figura 7 
A 2 saída 
mudança de setpoint, [  moceso | ——— 
pelo operador 


Identificação pelo teste da resposta em frequência; o pro. 
cesso é submetido a uma entrada harmônica (sinal senoidal) 
De acordo com as curvas de magnitude e fase, é possível iden. 
tificar as frequências de corte (avaliando-se a influência dos 
zeros e polos) e, consequentemente, a função de transferêncii 
estimada correspondente, conforme mostra a figura (2.3). 


Fegistrador 


Figura 2.3 - Modelagem pela resposta em frequência 
analisador de 
espectro 


BEE a [ o | 


IMagíjc)| IFase(jo)] 


Pt É prec 
OS € zeros e, eventualmente, o val” 


las! 


dereeses Oz» mm a 
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do atraso de transporte ou tempo morto. Todo o cálculo dos 
parâmetros é feito em um momento e o tempo de avaliação é 
diferente daquele em que se realiza o ensaio. Por isso se diz 
que a estimação é feita off-line. Para este tipo de abordagem, 
utilizam-se na prática modelos discretos para os processos. 
Isso se justifica pelo fato de que os algoritmos de identifica- 
ção trabalham com os valores de amostras dos sinais de en- 
trada e saída. 


Figura 2.4 - Modelagem off-line 


disco 


[0] 
jerador de entrada pr iara 7 
| gera “eme processo o DAS, 


* Identificação on-line: é um procedimento iterativo via com- 
putador. A identificação offline tem uma desvantagem de 
implementação. Muita quantidade de memória pode ser 
necessária para armazenar os valores digitalizados de todas as 
amostras de entrada e saída do sistema obtidas no ensaio. 
Neste sentido, muitas vezes precisa-se de métodos recursivos 
que utilizem pouca memória e que sejam capazes de atuali- 
zar a estimação dos parâmetros do modelo a cada período de 
amostragem. Conforme representado na figura (2.5), diversos 
algoritmos estão disponíveis na literatura para realizar a 
estimação on-line. O mais popular é o algoritmo dos Mínimos 
Quadrados Recursivo - MQR. Em muitas aplicações, as medidas 
do processo são obtidas sequencialmente (capturadas a cada 
período de amostragem) e processadas on-line em algoritmos 
de estimação recursivos. A aplicação em tempo real dos algo- 
ritmos de identificação é interessante para vários propósitos, 
entre os quais rastreamento de parâmetros variantes no tempo, 
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tro 
ico, fi em, con 
detecção, diagnóstico, filtrag! 

' i ificials- 
vo e redes neurais artifici: 


rrqura 2.5 - Modelagem 9º 


entrada 


processo 


2.2 QUALIDADE DO MODELO MATEMÁTICO 
ESTIMADO 


A identificação, de modo geral, consiste em três etapas: deter. 
minação da estrutura, estimação dos parâmetros e validação do modelo, 
A figura (2.6) ilustra as diferentes etapas da identificação de processos. 

À realização do experimento proporciona a aquisição das medida: 
do processo (dados experimentais). Assim, depois de selecionado o perio- 
do de amostragem para cada processo em particular, coleta-se um con- 
junto de medidas (visando-se à obtenção de um modelo discreto), sendo 

Z=(y,U) 
Y=ly(D/t=1,. 
Ea ro NH; Usfu()/t=1,.,N) 
º conjunto de val, E 
conjunto de medidas Faia do sinal de saída, y(t), Ué? 
experimentação e t especifica os j Sistema, u(t), N é o universo 
Ti (To 2T, 3T, instantes de amostragem múltiplos & 
modelo, adequada ao sistema Po a determinação da estrutura & 
los Candidatos, Pois é realizada a partir de Ls 
i em eciment 
a seleção do conjunto de tificado e nas Ec oa e pa 
tificação ado : Ssores ções necessárias 
E Ae seia do tipo para PATA à aplicação da técnica de 
"A estimação do models útiea Ou não paramétrica, re, 
Psiste na seleção de técnicas 
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estimação para a classe da estrutura especificada. No caso da estimação 

amétrica, objetiva-se aproximar os parâmetros desconhecidos dos 
polinômios do modelo (representação entrada/saída). A partir de Z", a 
validação do modelo matemático é efetuada pela comparação do modelo 
obtido com as medidas em um conjunto de dados que não foi utilizado 
na estimação do modelo matemático (generalização do modelo). 


Figura 2.6 - Etapas do procedimento de identificação 


medidas do 
processo 


determinação 
da estrutura 


! 


estimação do 
modelo 


validação do 
modelo 


A estimação de parâmetros é um procedimento numérico que 
determina os valores dos parâmetros do modelo (desconhecidos) e pode 
ser formulada como um problema de otimização no qual o melhor 
modelo é aquele que se ajusta às medidas para um dado critério. Uma 
vez obtida a magnitude dos parâmetros, objetivando qualificar o de- 
sempenho do modelo estimado, técnicas de validação de modelos são 
aplicadas, Entre as diversas técnicas de validação, pode-se utilizar a 
comparação das respostas do sistema real com a resposta obtida pelo 
modelo estimado (pode-se, também, aplicar técnicas no domínio da fre- 
quência ou métodos estocásticos). O modelo é adequado se o erro come- 
tido no ajuste está em valores preestabelecidos, ou ainda, para uma 

aplicação, se a resposta do modelo estimado reflete cor- 
Tetamente a resposta do sistema, não apresentando discrepância entre a 
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a modelo (não existe um ms 

da do sistema reale a saída estimada pelo pao Pr mé á 
Es ificaçã é universalmente ace) E riado), 
do de identificação que justar Os parâmetros do modelo estimado 


O procedimento de ai : teens 
ia entre a saída do sistema € a do T j Mínimo, d 
Es ren está mostrado na figura (2.7). 


acordo com um dado critério, 
27 - Diagrama de identificação de ajuste do modelo 
igura Diagrama de ide 


saída 


— sistema 
entrada 
Emas 
E 
= ermo 
— [os | 


algoritmo de 
ajuste dos 
parâmetros 


9 cálculo para ajuste do modelo é reali ntrada e saída são amostrados é 
cursivamente, respectivamente, izado não recursivamente ou ré 


2.3 IDENTIFICAÇÃ 
CÃO E 
dia COMPUTADOR = RIMENTAL 


na figura (2.8), aprese” 


de d, 
Macenadaao8 de entrada e saída do 5% 
PO computador ou em o! 
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b) As medidas são processadas para estimar os parâmetros do 
modelo sem restrição de tempo computacional (utilizam-se algo- 
ritmos não recursivos). 


Figura 2.8 = Identificação de processos off-line 


entrada 


disco 
computador 
CT 


modelo 


Figura 2,9 - Identificação de processos on-line 


entrada processo saida 


E computador 
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Por ou 
figura (2.9), apresenta 
a) Poucas medidas da experim 
e processadas; en 
ritmo de estimação recursivo é utilizado para aut 


b) Um algo: à do de amostragem; 


os parâmetros após cada perí 
o) A quantidade de cálculos (tempo de processamento do algo. 
ritmo) necessária para ajustar o modelo é função do período de 
amostragem. 


Uma grande variedade de métodos está disponível na literatura 
para identificar processos, tanto on-line quanto off-line. As técnicas de 
estimação de parâmetros de processos podem estar apoiadas em 
sistemas de conhecimento (inteligência computacional), em tarefas de 
supervisão (auxiliar o usuário na utilização e regulação dos algoritmos)e 
em modelagem. Alguns produtos de software disponíveis comercial- 
mente para modelagem de processos monovariáveis e multivariáveis 
são: ADAPTs, AIDA, CONTSID, IDENT, PITOPS, SIMNON, WINEACT, 


PROTUNER, INTUNE, DELTAV, LOOPPRO, P ol. 
MONITOR e PROCESSDOCTOR (VANDOREN, 2008) aa 


2.4 BACKGROUND E LITERATURA 


A tomada de decisões e 


do acesso à informação Fer resolução de Problemas são dependente* 


a 
uada sobre o 


ervados é denom ? 'elevante descrição do 


cessos) ea descrição E) Nominada identi 
Via, existem diversas fo; sistema Tesultante é ntificação de sistemas 


A segui es Ê nominada modelo. 
Seguir, mencionam-se algumas dei que é identificação de proceso! 
* estabelecidas na literatu? 


dh 


dm 
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* Identificação é a determinação, baseada nos sinais de entrada 
e saída, de um sistema em que, para uma classe especificada 
de sistemas, o sistema sob teste é equivalente (ZADEH, 1965); 


* Identificação de sistemas representa a “interface” entre mo- 
delos do mundo matemático e o mundo real (JOHANSSON, 
1993); 


* Aideia da identificação de sistemas é permitir a elaboração do 
modelo matemático de um sistema dinâmico, baseado em me- 
didas coletadas pelo ajuste de parâmetros e/ou do modelo, até 
que a saída do sistema coincida, “tão bem” quanto possível, 
com as amostras das saídas medidas (LJUNG, 1996). 


A procura por modelos relevantes pode também iniciar sob outro 
aspecto. Considerando-se que a estrutura, ou projeto, de um sistema é 
amplamente conhecida, frequentemente é possível produzir um diagra- 
ma de blocos ou esboçar o rascunho do sistema e seus componentes 
funcionais (modelos mecânicos e elétricos). Tal prática é denominada 
modelagem. A complexidade da modelagem requerida depende do pro- 
pósito de modelagem e da identificação do processo, ou seja, verificação 
de modelos teóricos, sintonia de parâmetros de um controlador, projeto 
assistido por computador de algoritmos de con-trole digitais, controle 
adaptativo do tipo self-tuning control ou auto-tuning control, monitoração 
dos parâmetros de processos, detecção de falhas, otimização de siste- 
mas, entre outros. 

Outras importantes informações sobre identificação são: a) o tipo 
de modelo do processo (linear/não linear, variante/invariante no tempo, 
Paramétrico/não paramétrico, contínuo/discreto, monovariável/multi- 
Variável, característica do ruído, atraso de transporte); b) a precisão re- 
querida do modelo (baixa, média, alta); c) o método de estimação do processo 
(offHline/on-line, malha aberta/malha fechada, resposta ao degrau, respos- 

ta em frequência, análise espectral e de Fourier, técnicas de correlação, 
de estimação de parâmetros, etc.). 

À identificação e a modelagem de sistemas têm diversas raízes, 
nas mais diversas áreas do conhecimento, seja em matemática, esta- 
tstica, ciência da computação, teoria de controle, análise de sistemas, 
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essencialmente começa com K. F. Gauss em 1809 e termina por vol. ce 

ta de 1960, quando modelos paramétricos explícitos começam a (é 
Scasionar interesse na comunidade de controle. Durante esse perio- 

do, todos os conceitos estatísticos essenciais usados em sistemas ae] 
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9 1970-1985 (consoli Posta ao imp úlão an 
: ): as li 


entre difarentes Mestras ia Prática em sistemas de po 
Ware para q iria etodo) ogias e dat er distinguidas em comi 
Moda Problemas o desenvolvimento de Sof” 
em identificaçã O; 
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d) 1985- ... (surgimento de novas ideias em raízes estatísticas): em 
meados dos anos 80, a visão estatística de sistemas de iden- 
tificação já está madura e sedimentada. As estruturas clássicas de 
estimação de parâmetros alcançam relativo sucesso e podem ser 
coerentemente transferidas para o mundo dos sistemas diná- 
micos. O desenvolvimento de ferramentas computacionais pode- 
rosas e o aparecimento maciço de pacotes de softwares comerciais 
também contribuíram para o aprimoramento da área de identifi- 
cação. Entre os tópicos com desenvolvimentos significativos, hoje 
em dia, que na maioria das vezes possuem pequena interação com 
métodos estatísticos, citam-se: métodos de subespaço para mode- 
los espaço de estados, identificação para controle, propriedades de 
rejeição e tratamento de ruído, modelos black-box não lineares 
(área de redes neurais artificiais, modelos fuzzy, transformada 
wavelet, etc.) e aplicações e análise de dados no domínio da fre- 
quência (LJUNG, 1999; MARCHI et al., 1999; AGUIRRE, 2007). 


2.5 APLICAÇÃO EM CONTROLE ADAPTATIVO 


Um exemplo de utilização da identificação de sistemas em en- 
genharia de controle de processos que tem sido enfatizado nos últimos 
anos é na implementação de controladores adaptativos e/ou preditivos. 
Justifica-se a utilização destes algoritmos de controle pela adequação 
para tratar processos variantes no tempo com parâmetros desco- 
nhecidos e, adicionalmente, em aplicações práticas devido à ineficiência 
do desempenho dinâmico de algoritmos de controle com ganhos fixos 
na presença de complexidades (dinâmica assimétrica, parâmetros va- 
nantes, não linearidades) no processo controlado (ROFFEL et al., 1989; 
MALIK et al, 1991; GESSING, 1996). 

O controlador adaptativo do tipo autoajustável (self-tuning) apre- 
Senta uma malha adicional que automaticamente modifica os parâmetros 
Para compensar as variações ocorridas no processo ou no meio ambien- 
te, trazendo benefícios significativos em relação a outros tipos de 
Sontroladores, pois as variações da dinâmica são acompanhadas pelo 
identificador e, consequentemente, pelo controlador a cada período de 
amostragem. Entre os controladores adaptativos, o regulador autoajus- 
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48 EE E 
; sto por K: J. Âstróm e B. Witten 
tável de variância mínima proposto por * ara 
ELE tido crescente utilização, implementação e varia O d 
em 1973, le de processos (COELHO et aj, or 
projeto no contexto de controle E t d + 19 1 
ASTRÓM; HAGGLUND, 2006). A partir desta estrutura de controje 
possível a elaboração de controladores adaptativos robustos, isto é 
controlador autoajustável de variância mínima generalizada (c A 
GAWTHROP, 1975) e, posteriormente, o controlador preditivo gener, 
lizado (CLARKE et al., 1987). 

Esta estratégia de controle está dividida em duas etapas; iden. 
tificação e controle, conforme ilustra a figura (2.10). A etapa de identi 
ficação consiste na obtenção dos parâmetros do modelo do si 

Pp: Sistema 
controlado, através de um estimador de parâmetros recursivo. A segund 
etapa, de controle, consiste no cálculo dos parâmetros do controlador 
pela eee de um critério de desempenho especificado (procedi. 
mento ótimo) ou pela técnica de alocação de polos (procedime: 

3 nto clás- 
sico) (SEBORG et al, 1986; BOBÁL et al., 2005). 


Figura 2.10 - Cont 
gura 2.10 - Controlador adaptativo do tipo autoajustável 


entrada 


| modela do 


E O controlado, ing 
insere um relé no Es a Sea Autossi 
ordem reduzida a Partir da Sa estimação End spin a: 
Posta em frequência 
. Os experimente” 


NS 
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com relé na malha de realimentação em controle adaptativo tornaram- 
se populares em equipamentos industriais a partir do trabalho de 
Ástrôm e Hagglund (1984). Este método é utilizado para determinar o 
ganho crítico e a frequência crítica e, por consequência, automatizar 
os métodos de projeto (sintonia) de controladores PID (Proporcional- 
Integral-Derivativo) proposto por Ziegler e Nichols (1942). A figura 
(2.11) ilustra o elemento não linear (relé) realimentando um sistema de 
controle para os propósitos da identificação paramétrica ou calibração 
do controlador (VISIOLI, 2006; YU, 2006). 


Figura 2.11 - Controlador adaptativo com autossintonia 


referência + y sistema | saída 
EE 


relé 


Observação: 


* Os controladores industriais com autossintonia podem elimi- 
nar a função do operador ou engenheiro especializado, únicos 
capazes de obter resultados satisfatórios na sintonia manual 
dos parâmetros de um controlador PID. Aplicações típicas de 
controladores com autossintonia têm sido utilizadas, prin- 
cipalmente, na busca por melhores resultados operacionais 
quando comparados aos apresentados pelo conjunto de pará- 
metros atuais ou padrões de fábrica. A autossintonia é apli- 
cada, muitas vezes, como um procedimento de inicialização 
de um controlador adaptativo, denominado pré-sintonia. É o 
caso de controladores embarcados com a técnica gain schedu- 
ling, onde em cada ponto de operação do processo é definido 
um conjunto de ganhos do controlador. Um equipamento 
industrial é o ECA-600 do fabricante ABB, que dispõe das téc- 
nicas de sintonia gain scheduling, auto-tuning e self-tuning. 

Nasi 
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1) Tustrar uma experimentaçã 
modelagem baseadas nas respostas à 


offline e on-line. 


o e descrever à aplicação das técnicas 
o degrau e frequência, identificação 


modelo estimado pelos testes da resposta a 


2) Avali alidade do E 
valiar a quali uencial (diagrama de Nyquist), considerando 


degrau de malha aberta e freg! 


processo real 


uti 1020s+4080 
245! + 2185 + 9565? +2456s + 4080 
modelo estimado 
E Ee E 7.34 
s +1.68s+7.34 


3) Avaliar a qualidade dos modelos estimados em relação aos testes das 


respostas ao degrau de malha aberta e fr jal (diagrama de Nyquist) 
equencial : 
para o processo real e os seguintes modelos poda a 


processo real 
Dot e sen 
(0.62s + 1(17.25+ 1) 
modelos estimados 
[em (s) = (VBA ess liso 
1.98 + Da? 
À cu É 
Cls)= AR 
4. 
ua S+1(167571) 


tamento temporal de mo Modelos 
tário utilizando um conaha fechada parádos com base no compor 
sintonias,istoé, — “Strolador py (pol2 Uma entrada degrau 
' 1 Proporcional. tora) ué 
HIntegr. 


A caractt 
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4) Uma malha de controle com controlador PI, uma válvula não linear e 
um processo de terceira ordem estão ilustrados na figura (2.12). 


Figura 2.12 - Malha de controle analógica com controlador, planta e válvula 


E U MV 


A característica da válvula está representada por 
v=f(u)=u” ; u>0 


Verificar por simulação (Simulink) a necessidade da ressintonização do 
controlador PI ajustado com K. = 0.15 e T; = 1, assumindo os seguintes 
níveis de operação: 

a) Mudança de referência de 0.2 para 0.3; 

b) Mudança de referência de 1 para 1.1; 

€) Mudança de referência de 5 para 5.2. 


À sintonia PI deve garantir, nos três níveis operacionais, uma dinâmica 
de malha fechada sem sobressinal e com tempo de resposta de aproxi- 
madamente 10 s. Idealizando-se uma concepção de controle do tipo PI 
ganho escalonado (gain scheduling) e pela sintonia tentativa e erro, quais 
9s ganhos nas três faixas de operação? 


9) O atraso de transporte pode ser real (transporte de material ou ener- 
ou o efeito causado pelo acúmulo de várias constantes de tempo. 
aproximações polinomiais podem ser utilizadas para representar 
9 atraso, ou o atraso é usado para representar a resposta de um sistema 
de alta ordem. 
Bites 
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52 e E 
a) Obter K, re 0 por simulação, onde as pi & 
planta, G,(s), e do modelo, G,(8), têm ae 
' ) j 
G=no é 60=qa+T 


POE (ora)! 


b) Encontrar 7 por simulação ou epa oii as fu. 
ções de transferência da planta e do modelo são exp po 


eos x 1 
Go(9)= 3541 * Go(8)= (3s+1)(ts+1) 
1.66e 


j : : G,()=>——. Estimar 
6) Seja a seguinte planta controlada: G,(s) (1,0754517 por 


simulação um modelo de primeira ordem com atraso de transporte d: 
A Ke 
forma G, (s)= e Observando a rapidez da resposta e a energia do 


U(s) 
G(s)=— 
“= 26) direta 
de sintonia 7=4,8, do oi 
1+% 
K.= 
K(t. +04) , T=1+0% F brasão 
2 


conto 
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MODELOS 

DE PROCESSOS DE 
ORDEM REDUZIDA 
E COMPLEXOS 


3.1 INTRODUÇÃO 


Vários processos industriais complexos são aproximados de for- 
ma apropriada por funções de transferência de baixa ordem, ou seja, 
modelos matemáticos de primeira ou segunda ordem. Os sinais po 
dem ser representados utilizando-se variáveis contínuas ou discretas e 
as funções de transferência, através das transformadas de Laplace e z, 
respectivamente (SEBORG et al., 2004; ÁSTRÔM; HÁGGLUND, 2006; 
OGATA, 2010). 
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«1 
3.2 MODELOS CONTÍNUOS DE ORDEM REDUzinA 


A re) is ust um modelo matemá,; 
presentação dinâmica mais US al de E e ático 

uma planta industrial é à função de tr: e eir 

para s 


ordem com atraso de transporte, isto é, 
Os 
Y(s) Ke 
Gal) = O) teia eu 
pe U(s) ts+1 


onde s é o operador Laplace (s=d/dt), K, é o ganho, O é o diacto, de trans. 
porte contínuo, 7 é a constante de tempo e Y(s) e U(s) são as transfor. 
madas de Laplace da saída e da entrada, respectivamente. 

A seguir, discute-se como os parâmetros (K,, O, 7), que caracte. 
rizam o comportamento dinâmico de um sistema de primeira ordem, 
podem ser identificados com base na curva (resposta) de reação do 
processo. 


Ganho de regime 
O ganho de um sistema é definido por 
K AD y(o)-y(0) = Yt)-y(0) 


P 


* Áu(t) é à variação da entrada e té 


determinação do porão, O exemplo (3.1) mostra o procedimento de 
j Stro K, para um processo de nível. 


tável, conforme a figura GI) on de nível de líquido em regime es 


cálculo do ganho conside: tram os comi ' 
po; Jeo 
Brau) na vazão de entrad pr Um aumento dios nie a 


(através da um a 
S da abertura da Válvula de saí “pis degrau na vazão de saída 
É Ctivamente, 


<R 


as 


[E 
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Figura 3.1 - Sistema de nível de Líquido 


entrada Mi 


nível 
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IpENTIFICAR ES = Ea 


co a 


Observações: os das figuras (3.2) e (3.3), Observa, 


qultad 
. ordo com os rest iza uma mudani 
De e e E quando se realiza u ça No pro. 
seq 


cesso; 


* Oganho de um pi 
ser avaliadas cuidadosamente par: 


rocesso tem unidades específicas que devem 
a cada aplicação particular, 


Constante de tempo . 4 
A constante de tempo para sistemas dinâmicos (caracterizada em 


uma função matemática de primeira ordem) é definida pelo 
tempo requerido, medido a partir do ponto onde o sinal de saída 
começa a mudar, para que a saída do processo de primeira ordem 


atinja 63.2% do valor da variação total, depois que ocorreu uma 
mudança na entrada. (MOLLENKAME, 1988). 


* Por que 63.2%? 


Para um processo de Primeira ordem e sem atras: 
o de transporte, à 
equação diferencial é ga 


dy(t) 
T E +y(b= K,u(t) (3.3) 
onde y(t) é a saída, u(t) é 
fenda Ee ai K, € 0 ganho e ré a constante de 


m degrau de amplitude M é dada por 

YO) =K Ma -e%) : 
Portanto, observando-s 
Para valores múltiplos da constant dotes: à Saída pode ser avaliada 
º* yD= 0.632K,M » isto é, 

* yQ)= 0.865K,M Tresponde à Constante de tempo) 

* y(37)= 0.950 

* y(47) = 0,989 


* Y6)=0.993k Mm Corresponde 2º tempo de acomodação) 


t>0 (34) 
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Atraso de transporte (tempo morto) 

O tempo morto na dinâmica do processo normalmente é o tempo 
que leva para o material se mover ou ser transportado de um ponto a 
outro. Esta é a razão pela qual o termo equivalente “atraso de trans- 
porte” sempre é utilizado para descrever tempo morto. 

O atraso de transporte pode ser definido como 


o tempo decorrido após a ocorrência de uma perturbação no 
processo até que seja notada uma mudança na saída do mesmo. 
(MOLLENKAMP, 1988). 


O exemplo (3.2) descreve a presença do atraso de transporte num 
processo de mistura de fluidos. 


Exemplo 3.2 — Dois fluidos, quente e frio, são misturados num sistema 
de bombeamento conforme ilustra a figura (3.4). A temperatura das 
vazões combinadas é medida em um ponto a jusante do ponto de 
mistura. Qual é a resposta da temperatura medida após um aumento em 
degrau na vazão do fluido frio? 


Figura 3.4 - Mistura de dois fluidos, quente e frio 


quente 


sensor 


> 


vazão de fluido frio 


temperatura no ponto de mistura 


temperatura medida 
no sensor 
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à 
bem representados por uma funç 


E isto é, 
de acordo com as equações (3.5) ou (3.6), 18 
-0s 
Kk,e 
E +1) (35 
Sora rt +) ) 
k,e de 
Goals) = E R 
re RAE ea! 


onde 1; e ” são as constantes de tempo, K, é o ganho, O é o atraso de trans 
Porte contínuo, (v, é a frequência natural e S éo fator de amortecimento. 


* Sec ; E: sistema sobreamortecido (raízes reais e desiguais, 
nn); 
* Seç=1 sistema criticamente am rte 
asas Ortecido (raízes reais e iguais, (| 
os 3 
É <T:sistema subamortecido (raízes complexas conjugadas). 
Observações: 
Seja a fu; ET 
poli eo do dopalrênci (9 UO=BO/a(. pipe 
Sistema, en, uanto que as RA 4 pan denominadas aaa É 
(), são denomi, i 


Polinômio do numerado!; 
ma; 


istema de terça; 
Buint; er, 
e função de transferência,” Ordem Caracterizado pela se 


rp 
Cos 


n(ys+1) 
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Se a condição |1/)|> 10/60, é assegurada, então a resposta do 
sistema de terceira ordem pode ser aproximada pelas raízes 
dominantes do sistema de segunda ordem (DORE; BISHOP, 


1995); 


* A dominância entre polos contínuos é dada pela relação 10: se 
ps é o polo dominante de p;, então |p,|>10)p;|. A relação polo 
real dominante no caso discreto é dada por |p,|>19/p, 


* A presença de um polo próximo da origem no plano-s e no 
Plano-z proporciona respostas temporais lenta e rápida, res- 
pectivamente; 


* Modelos matemáticos de planta são adequados para o projeto 
de controladores, por exemplo, o controlador por alocação de 
polos. Uma prática comum é sintonizar os controladores PI e 
PID com modelos de ordem reduzida, nestes casos, de primei- 
ra e de segunda ordem. 


3.3 MODELOS DISCRETOS DE ORDEM 
REDUZIDA 


Embora muitos processos sejam contínuos por natureza, os 
modernos sistemas de controle utilizados para controlar processos 
co Se em computadores digitais e aplicam algoritmos de controle 
digital (JACQUOT, 1995). O sinal de controle é implementado em in- 
E s de tempo com o intervalo denominado período de amostragem, 
5: Se a técnica de Projeto de controle por computador baseia-se no 
ra matemático do processo controlado, então um modelo discre- 
Planta é necessário para Projetar e calcular a saída do sinal de 
Fig As diferentes formas de representar uma malha de controle 
na figura (3.5). 
int 


ae 
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YiS) +, a “649 
Fis) j 


(a) sistema de controle « 


YAS) Gs) RA sa 
(s) + 
CE | Eq 


(b) sistema de controle digital (por computador) 


" Ye) Gal 
CIO| “8 |—| go) 


(c) sistema de controle discreto (amostrado) 


ontínuo (analógico) 


> 
É 
E 
[=] 
[4 
Ê 
Ê 
Ê 
E 
e 
Ê 
E 
E 
Ê 
E 
= 
e 
aê 


conjunto de amostras da saída a contínua do processo, * U” 


está disponível para identi Processo em períodos de amos fm, 
(3.6) j e controle i 
), OU seja, + Conforme ilustra à 
E 
Figura 3.6 - Segurador de ordem ze; q 


Bric 
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saída contínuo e pulsado (discreto) e U(s) e U'(s) são os sinais de read 
contínuo e pulsado (discreto). 
A função de transferência amostrada da planta é calculada por 


1109) Go(8) 
votos ont am am 
e sendo Gy(s) da forma 
=T, 
is 6) (3.8) 
s 


A função de transferência discreta correspondente à equação 
Gp(s), equação (3.1), é 


zºb, E po aD, 


z+a, l+a,z” 


(3.9) 


onde a; e bo são os coeficientes do modelo matemático discreto, z é o 
operador deslocamento, z”y(t) = y(t-n), d é o atraso de transporte discre- 
to e satisfaz a seguinte relação: 


dT, <O<(d+I)T, (3.10) 


O parâmetro a; é função de re T,, enquanto que o parâmetro by é 


função de K,, 1 e T;. 
Por outro lado, a função de transferência discreta correspondente 
à equação G,:(s), equação (3.5), é 


Gts 2“(byz+b,) “Db +bz!) 
E = 


- (211) 
2 +azra, lraz'+az 
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Com o auxílio da tabela de transformada z ou com um software 
aplicativo em controle (Matlab), é possível obter, a partir da função de 
transferência de processos contínuos de ordem superior, a equivalente 
função de transferência discreta. 


Observação: 
= As funções de transferências discretas de processos contínuos 
de primeira e de segunda ordem, incluindo o segurador de ordem 
zero, são dadas pelas equações tabela (3.1) (ROFFEL et al, 1989). 


Exemplo 3.3 — Considere o processo de primeira ordem representado por 


K 
G;()=—E 
ps) ts+1 


onde Ky = 1, T= 7.55, T;= 4s. Obter a função de transferência discreta 
sem e com o segurador de ordem zero. 


Ko K 
ts Se E 
G,(s) E Samos Ko pl 


a) Caso sem o segurador de ordem zero (a estrutura do modelo é não causal) 


Eee ba K 
nega RO 4 


a,=-e** =-0.5866 


( =0133 ; 


b) Caso com o segurador de ordem zero (a estrutura do modelo é causal) 


boz”* 


G'(z)= 
49) l+a,z! 


RA 
bo=(1-e"%) "A =0,4134 ; a, =-e"! =-0.5866 


Es Bintore 
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: tros). 


estrutura e nos parâme 


Observação: ta do sistema discreto de Primeira 


: é q 
“A epa Si rice Re comportamento r a 
inss quando o polo real -a: dera por pa, 

é i iva ao círculo unitár k) 

plo, do ponto 1 até a origem relativ: Unitário er fe 


plano complexo z (FADALI; VISIOLI, 2009). 


3.4 MODELOS DE SISTEMAS DISCRETOS 
COMPLEXOS 
Se um computador (implementação digital) é utilizado pn 


controlar um processo, o sistema de controle de malha fechada pode se 
representado pela configuração apresentada na figura (3.7). 


Figura 3:7 - Sistema de controle por computador 


Conversor digital/analógico 
o 


modelo enj 
daz, é obtido por 


trada/ 
pg Processo, utilizando a transform” 


Y(z)= So(2)U(z) qui. 


à função 
dem zero, po nferên 


e 
CM 27, isto é 


onde G(z) é 
rador de or; 


Polinômios cia 


de ser repre 


sa Processo, incluindo 0 e 
“entada por uma razão de 


pa 
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Da BS) Yi) 
G,(z)=2 Ag) U(z) (3.13) 
ou 
A(zy(t)=2 “Blz )u(t) (3.14) 


A equação (3.14) é uma representação usual na literatura de 
Identificação de Sistemas Dinâmicos e Controle Digital Direto. Os polinô- 
mios A(z7) e B(z7), relacionados aos polos e zeros de malha aberta do 
processo, respectivamente, são 


A(z))= l+a;z” +az? desashá ge 
Blz )=bo+b;z" +b,z? +... +bz 


A correspondente equação a diferenças da planta é 
y(O+ay(t-1)+...+a y(t—na)=bou(t—d)+ 
bju(t-d-1)+...+b,u(t-d—nb) 
onde y(t) e u(t) são os sinais de saída e controle, respectivamente, dis- 
poníveis nos instantes de amostragem t = nT;, onden € Zº. 


Exemplo 3.4 — Obter a equação a diferenças equivalente e avaliar as ca- 
Tacterísticas de malha aberta do modelo discreto relativo ao processo de 
segunda ordem representado por 
gar ca 

(3s+1)(5s+1) 


Considerando o segurador de ordem zero e T, = 1 s, então com o 
auxílio do Matlab tem-se 


G,(9) 


y(=1.5353y(t-1)-0.5866y(t — 2) + 
0.0279u(t — 1) +0.0234u(t — 2) 
De acordo com a equação (3.14), observa-se que 
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3z 
aqu 1)=1-1.5358 


Blz 


Assim, como as raízes rela! 
interior itário, no Pp! 
interior do círculo uni fo 
de fase mínima. Os resultados do exemp! 


cutando-se o programa em Matlal 


Tabela 3.2 - Programa em Matlab da conversão contínua para discreta 
abela 3.2 - 


IDENTIFICAÇÃO DE SISTEMAS DINAmiça 
24 
=. 1, 0.5866z * ; 


1)=0,02794 0.0234z 


tivas aos polos e zeros posicionam. 
lano complexo z, O sistema é es 


b da tabela (3.2). 


s UM 
Y 


1. d=1 


(3.4) podem ser obtidos e, 


rte uma função de t 
(s)/den(s) - 
onv([3 1],[5 
2dm (num, den, Ts, 'zoh') 
do sistema discreto!) 
oots (numz) 


olos do sistema discreto!) 
roots (denz) 
(numz, denz) 
title('resposta a se 
XHabel ('satda"),xlabel (1 


amostras ') 


ransferência contínua para discrata 
> Gp(z)=numz (z) /denz (z) 
1]);Ts=l; 


quência unitária” ) 


&. 


% Modelo discreto 


% Zeros e polos 


% Resposta 


PESTE E 


DISCRETOS 

Aperturbaçã 

ils Pode ser definida como 

que tende a af, 

ef etar ad, 

As co nd EN 
Seguinte figos vinte pi Composto do processo é dá E % 

te figura (38) “da, de processo e 


SPSS S 
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cxpituLo3 Model 


Figura 3.8 - Representação do processo com perturbação de carga 


perturbação 
processo v(t) = cte 
uít) Es B(z!) y(t) 
A(z!) 


cuja representação matemática da planta é 


B(z!) 
A(z!) 


v(D= u(t-d)+v(t) ; vt)=constante (3.15) 


Neste tipo de representação da figura (3.8), é habitual distinguir, 
além da perturbação de carga, diferentes tipos de perturbações ou 
imprecisões que afetam a modelagem e estabilidade do processo con- 
trolado, tais como: perturbação estocástica, efeito não linear, dinâmica 
não modelada, erro de medição, variação paramétrica ou imprecisão na 
estrutura do modelo da planta (modificam a saída e não estão em G;(2)). 
Quando presente em um sistema de controle, o offset é eliminado pela 
ação integral do controle, mas alguns esquemas de controle adaptativo 
estimam o valor de v(t) (amplitude da perturbação) para compensar seu 
efeito e melhorar a dinâmica de malha fechada da planta. Entretanto, as 
diferentes formas que o sinal v(t) pode assumir são: pulso, degrau, 
*ampa, senoidal ou ruído aleatório. O projeto dos algoritmos de con- 

a trole ótimo depende da natureza da perturbação. Na teoria de controle 
p “stocástico, os controladores são desenvolvidos pela suposição de que as 
qêr, são estocásticas por natureza. A forma importante de 
(aê? Poti aleatória é aquela associada com mudanças não previsíveis 
Sistema. Estas perturbações são agregadas ao processo e modeladas a 
Partir de ma sequência ruído branco (ruído Gaussiano), e(t), com média 

5 Sé, conforme ilustra a figura (3.9). 
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TIFICAG: 
IDEN 


â nitária 
la e variância u 
ula 
= dia n 


mé 
com 
neo 

4,9 Ruldo bra 

Figura 3 


*=rand(N,1); 
Ye2ta-1; 


SR Ág) % Mégia 
eva * Média n 
Pan (e) var (e); 
Plot (e) 
Ylabel (tg 


(3 xtaber E 
Figure (2) 


(ac e lagolo 
*label (tauto 


Mostragi) 
CCOrr (g, 9, N 
Corre] a, re 


er 
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ca 


Observações: 

» A determinação da média, variância e valor médio quadrático 
de uma variável aleatória uniformemente distribuída x que 
tem densidade de probabilidade 


Ew)= ae! paraxe[a,b] 


0, paraxe[a,b) 


é calculada pelas seguintes relações (BARBETTA et al., 2004): 
b b 
E(xy)=%= [x£.G)dx ; V(x)=0; = [xP E dx ; 
E | a 
EG)= [xºE,G)dx = V() - [EF 


Estes resultados podem ser aplicados na tabela (3.3) para a ge- 
ração das variáveis aleatórias x, y, obtidas a partir do coman- 
op! do rand do Matlab; 


= Um processo estocástico v(t) é denominado ruído branco se 

suas realizações no tempo estão descorrelacionadas, isto é, a 

| função de autocorrelação, E(V(t) v(t+0)), tem valor não nulo 
E somente para ó = 0. Observe que isso implica Efv(t))=0, vt. O 
* processo é denominado ruído branco Gaussiano se, adicio- 
malmente, em cada instante de tempo a função densidade de 
é Gaussiana. Caso contrário, o ruído é denomi- 
colorido (ASTRÔM, 1970; WELLSTEAD; ZARROP, 1991). 


admitindo-se que v(t) é da forma 


cz") 
A(z!) 
ser reescrita como 


v(t)= 


e(t) 
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à Ce) er) 
y()= e u(t-D aço?) Gay 


E: a(z Dy(t)=2 dp(z Hut) + Clz Dele) 


blocos representando o sistema de controle de ma 


i d ' 
O diagrama de esso e perturbação) está ilustrado na figura (ag, 


lha aberta completo (proc 


Figura 3.10 - Representação do processo com perturbação estocástica 


Tt “as 


e(t) c(z?) 


5| Do! 
A(z?) 
v(t) 
u(t) 
=; 


ziBtE) YO 
| Ca rd 


traso de tran últiplo 
Sporte que é um mi ko 
Sões estocásticas. PMostragem, quando sujeito a 


Es 


O modelo 
da equação (3 1€ 


; 
cs 
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cal - 


A(z Dylt) =2 “Blz Du(t) + ds e(t) (317) 


ou 


Atz Nay(t =2 “Blz ")Au(t) +Clz "e(t) 


A concepção da representação do processo pela equação (3.17) é 
importante para assegurar uma parte integral na lei de controle digital e 
proporcionar erro nulo em regime entre a saída e a referência para qual- 
quer setpoint e perturbação de carga com comportamento em degrau. 

De acordo com a equação (3.17), diversos tipos de modelos mate- 
máticos lineares discretos, mostrados na tabela (3.4), podem ser utiliza- 
dos em identificação e controle digital (variância mínima, preditivo). 


Tabela 3.4 - Família de modelos para representar processos dinâmicos 


2 
Fa 
Considere as seguintes notações da tabela (3.4): 
5 = Modelo FIR (Finite Impulse Response); 
' 
e” AR (Auto-Regressive); 


- Modelo Ê 
É MA (Moving Average); 


ARMA (Auto-Regressive Moving Average); 

CAR (Controlled Auto-Regressive); 

CARMA (Controlled Auto-Regressive Moving Average); 

CARI (Controlled Auto-Regressive Integrated); 

; (Auto-Regressive Integrated Moving Average); 
CARIMA (Controlled Auto-Regressive Integrated Moving Average). 
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TIFICAÇÃO DE 515 "STO DINÂMICOS | 
DENTE eg 

it 


Observações: i s e Controle Adaptat; 
“ncação de Sistema: vo ( 

E sido pe implementação dos controladores p 
ted os coeficientes dos polinômios A(z), Bley, 
ca) são desconhecidos, podem ser ep 2 Vad 
no tempo e devem ser estimados por técnicas de identi 
em uma particular aplicação experimental (MALIK et al, 199. 
SHAH et al., 1991; LEVINE, 1996; COELHO et al., 1999); 

* Aliteratura de Identificação de Sistemas emprega no lugar dos 
modelos CAR, CARMA, CARI, CARIMA as notações ARX, ARMAR, 
e onde X denota exogenous ou external input (LJUNG, 


" ne síntese, os modelos lineares discretos baseiam-se no mo- 
o geral da tabela (3.5), com as suas respectivas definições. 


Tabela 3.5 — ifi 
a 3.5 - Modelos identificados Por função de transferência 


Modo Ay) BE o CA 
| Rem t-+ CED ao 


D=D(z4)=1 
F(zt) =D(21)=1 
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caPÍ 


3.6 PROBLEMAS 
1) O sistema de controle da figura está sujeito a uma perturbação a 
partir de uma fonte e(t). Determinar a função de transferência discreta 


do sistema de acordo com a estrutura CARMA. Admitir um período de 
amostragem de 1 se segurador de ordem zero - ZOH (Zero-Order Hold). 


A(z)y(t)=B(z Dult-1)+C(z7)e(t) 


Ts 
E E O 
s + 
E 


2) Considere o sistema de controle digital ilustrado na figura. 


Fo [E º 


A função de transferência de malha fechada é dada por 


Y(z) 27(1.125+0.4552") 
% Y(z) (1+0.7577! +0.4552 *) 


3) Plotar Os polos e zeros no plano complexo z e avaliar a estabilidade 
fechada; 


a função de transferência de malha aberta (Gy(2)"? É o proces 
e de fase não mínima? Por quê? 


a diferenças do processo (polinômios A(z), B(7") e 
? 
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SU | 
Ss ty, RA 
Rn to ra 
E imação numérica da deriva ; º 
3) A partir da aproximaç a yin o] am) A a 
ávtaro) gl, a ç 
de segunda ordem ] 
ão a diferenças É 
em E rsireiranmeid B(2”') e atraso de transporte discreto) gi 
direto 
R L 
RA! E 
po 
4) Seja o processo representado pela seguinte função de t frtuaa 
0 
e =D K,(1+Ts) axnldord 


Us) (1+TIA+Ts) 


onde K = 1, T1= 105, T, = 
deve ser controlado por co) 


b) As alterações na função de tr 
COPO Apresentar um atraso de tra, 

9 A representação do 
(os polinômios Ag? De 


), 
5) Para o e 
relacionando ai dê Controle da 


u(t) | 


Ny 
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6) Considere o projeto dos controladores P, PI e PID para manter um 
nível desejado Hot) no segundo tanque de um sistema de controle de 
nível de líquido em dois tanques acoplados, conforme apresentado na 


figura (3.11). 


Figura 3.11 - Sistema de controle de nível com dois tanques acoplados 


tod sy— A 
j 


—— Hg 


O controlador deve comparar o nível atual (saída medida) H:(t) = y(t) com 
o nível desejado H»i(t) = y/t) e utilizar a diferença, ou erro, e(t) = Host) — 
Hut) para automaticamente alterar a taxa de fluxo de entrada Qt) = u(t) 
para o sistema por meio do controle da válvula V;. O correspondente 
diagrama do sistema está ilustrado na figura (3.12). 


Figura 3.12 - Diagrama do PID paralelo para o sistema de nível 


Controlador PID 


o: Els) Y 
2:32 [5 -08[empe 


tico que representa a dinâmica do processo linear e 
o tempo de segunda ordem de malha aberta é 
s- 
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IDENTIFICAÇÃO == SU, 
N 


É ng 
GA 2 +68+5 


dequado período de amostragem e obter sá 
” cjonar um dé racasto; 
s eras matemática discreta CAR do p » 


de malha fechada 
dente representação 
b) eira iii controlador PID e processo) a partir Pp 


guração contínua da figura (3.12); 

& Para o sistema de controle discreto do item (b), admitir uma 
unitária (referência) com as diferentes sintonias PID da tabela (36), 
avaliar o comportamento dinâmico do processo de malha fechada A 
mular com o Simulink) em função das especificações de desempenhe 
da tabela (3,7), Comentar os resultados da simulação digital. 


Tabela 3.6 - Sintonia dos controladores P, PI e PID 


7) Avaliar Por sim 


discretos das no 
sentada por CQuações 6 gra o adalos lin 

OMS (0) cara à 
S(9)=.1 


(+ 


“ Digitalizado com CamScanner 


CAPÍTULO 4 


MÉTODOS 
CLÁSSICOS PARA 
MODELAGEM 

DE PROCESSOS 


- 44 INTRODUÇÃO 


Na análise e projeto de sistemas de controle, deve-se adotar uma 
base de comparação entre os vários sistemas avaliados. Esta base pode 
ass especificando-se os sinais particulares de entrada e compa- 
as respostas dos sistemas visando, por exemplo, o menor 
de resposta. Os sinais de teste de entrada comumente utilizados 

degrau, rampa e senoidal. A figura (4.1) ilustra o procedi- 
avaliação da dinâmica de um processo de malha aberta perante 
teste, 
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E SISTEMA 
IDENTIFICAÇÃO DES 


78 SE. erante sinais de teste 
1 - Avaliação de processos peré 
il lt ql, saída 
entrada ooesso RR > 
SRS EN dinâmica do 
sinal de processo 
teste 


Para a identificação das características o um Processo 
sob avaliação, é possível utilizar Gs pe dg * raia [o Chaves, 
mento abrupto através de um acréscimo ou decr: e E Magnitude 
do degrau pode ser estabelecido pela variação da tensão (ou Corrente) 
pela abertura ou fechamento de uma válvula de entrada. Admitindo uma 
malha de controle, a resposta ao degrau pode ser obtida de acordo com 
Os seguintes passos (MOLLENKAMP, 1988; OGATA, 2010): 


a) Ajustar o regulador para o modo manual; 


b) Modificar a magnitude da variável de controle (acréscimo ou de- 
créscimo); 


O Registrar (plotar) a variável de saída do Processo. 


Alguns exemplos tí 
à resposta ao degrau de malha aberta (curva ER de 
cessos industriais reais) são mostrado: Teação da planta de p 

Às características de um sistema i AE dc ; 
trica, fase mínima, integrante) pod (estabilidade, dinâmica assimé 


i u é um : 

caracterizar, de fi Jimi Procedimento adequado par? 
: na pio ie do processo por meio à 

“tal, 2004; ASTRÓM; HÁCGLUND, 2000 JOSEPH, 2002; SEBOR 


Observação; 


de " linear Pela determinação da sº” 
ei tes amplitudes no sinal de e 
forma (curva) da resposte SIçÃO): “um sistema é linear S€! 
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RO 
capfrutos Métodos clássicos para modelagem de processos E) 
Figura 4.2 - Respostas típicas ao degrau de malha aberta 
saida 02 saida 
04 A E 
02 0,1 
: tempo o tempo 
o 20 40 60 0 20 40 80 
saida 
to saida 4 
D 
2 
5 
o 
2 tempo 
Y . 
o 20 40 E) 0 20 0 80 
saida 012 saida 
o2 0.06 
o 
-0,06 
% 20 40 E) º 2 bi E 
Macaca 
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80 Ee — 


é imento ade: 
Assim, a resposta ao degrau é um copiada diado par 
caracterizar, de forma preliminar, a dinâmica do P Por meio 


simples interpretação gráfica (BROSILOW; JOSEPH, 2002; SEBORG d 
al., 2004; ASTRÓM; HAGGLUND, 2006). 


Observação: | 

* Um sistema é conhecido por linear pela determinação da sua 
resposta ao degrau para diferentes amplitudes no sinal de en 
trada (teorema da superposição): “um sistema é linear se a 
forma (curva) da resposta ao degrau não depende da ampli 
tude de entrada”. 


A seguir, avaliam-se alguns procedimentos para a determinação 
dos parâmetros dos modelos matemáticos de um processo a partir dos 
compotamentos transitório e de regime permanente (métodos deter- 
minísticos, ou seja, métodos clássicos da engenharia de controle de pro- 
cessos e obtidos por inspeção da resposta da planta). Diferentes métodos 


Um modelo paramétri é 
dó dita dum pr cm 
transferência: tizado pela seguinte função 


Ta E —— (43) 


três p; , 
Ra Banho estático K,, constante Z 


ansporte 6 
Laplk ída à Sendo 
e da saída e entrada, Tespectivamer os € U(s) as transformadas 
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neirae 


sposta 
uência 


edidas 


TA 


mer” 


« Métodi 


Observações: 
* Asaída de diversos Processos químicos exibe Uma resposta ao 


degrau que pode ser aproximada pela equação (4.1), também 
denominada FOPDT - First-Order Plus Dead-Time; 


* O modelo discreto da equação (4.1) é expresso na forma 


z“B(z1) Y(z) 
Gln==—"== 2240) 
ni Az) UM 


y(b=a,y(t-1)+ bou(t-d-1) 


A(z")=1-a,27;a, =exp(-T /1); 
B(z!)= b, =K(1-a,); d=0/T,; 
* A qualidade dos parâmetros estimados da equação (4.1) está 


associada às medidas da planta cuja aceitação depende se a re- 
lação sinal/ruído > 5; 


* O modelo matemático da equação (4.1) pode ser utilizado na 
sintonia de diversos controladores: PID, IMC - Internal Model 
Control, alocação de polos, Dahlin, tempo mínimo, entre outros. 


*Arresposta ao degrau unitário (u(t) = 1,t >0) do modelo mate- 
mático representado pela equação (4.1) está caracterizada por 


YO=K,(1-e (4.2) 


aj base ma resposta temporal, é possível mensurar pontos para à 
“Plicação de diferentes métodos de estimação. Na literatura, existo uma 


de técnicas baseadas na resposta do processo ao degrau para 


tificação de K,, Te 0, Entre os disponíveis na literatura, avaliam-se os 
ke puintes métodos, apresentados por Ziegler/Nichols (1942), Sundaresan/ 


(1977), Nishikawa (1984), Smith (1985) e Hágglund (1991) 


OLLENKAMP, 1988; DORF; BISHOP, 1995; SEBORG etal. 2004). 
ma os métodos de Ziegler Nichols (ZN) e Hágglund (HAG), os parâmetros 


São calculados conforme ilustra a figura (4.3). A reta traçada cor- 
à tangente no ponto de máxima inclinação da curva de reação 
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E reação no fo 

ta da saída tem uma curva A de Saio emy) 
ip étodos, é calculado pelo : eMpo enç 
O atraso à, nos dois m "da do processo e o instante t, em 


do degrau na entrada e 
pi (t) = (0), No método de ZN, a constante de 


te 


a aplicação 


reta tangente toca à reta Pi 
E determinada pelo intervalo de tempo entre t; e 0 instanp, ; 


tempo 7 Je x 
em que a reta tangente toca a reta vo) =Yo ErquEnTO, que, nO Método 
HAG, ré calculada pelo intervalo entre ty e o instante t> em que à cur 
de resposta alcança o valor y(t) = (0) + 0.63yr. O parâmetro K, é cala 
lado pela equação (4.3), isto é, 

K =  HD-y(0) 

* Au u(o)-u(0) (43 
Figura 4 3 - Métodos de ZNe HAG para à modelagem de processos de Primeira ordem 


ví) 


ví) 
"Has 
3 
ES denis 
Método (il 
Yados og j de 5, 
resposta o tântes ith, 
Pelos pomposo po +, é Stva de 
210 2mgye “º Corregp o SÃO da planta são obst” 
“Na ey) Potdentes às passagens 
2), respectivamente 
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Métodos clássicos para modelagem de processo 


onforme apresentado na figura (4.4). Os parâmetros são obtidos pelas 
[e 


seguintes relações: = a 
"hu (4.4) 

t=1.5(t, —t;) EA 

0=t,-T=1.5t, —0.5t, ui 


Figura 4.4 — Método de Smith para à modelagem de processos de primeira ordem 


aii método de Sundaresan/Krishnaswamy (de forma semelhante ao 
qse de Smith) evita a utilização do ponto de máxima inclinação para 
A dnação dos parâmetros do modelo paramétrico. Dois instantes de 
Pri t; e to devem ser estimados a partir da curva de resposta ao 
Paes correspondentes aos tempos de resposta 35.3% e 85.3%, res- 

cp como ilustrado na figura (4.5). Os parâmetros são cal- 
ay (47) 
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; . 
7=0.67(t,— ti) 9 
0=1.8t, -0.29t, 9 
Figura 4.5 - Método de Sundaresan/Krishnaswamy para a modelagem 
de processos de primeira ordem 
vít) 
853% 
35.3% 
(O) 
ai Proposto por Nishikawa, É 
cálculo das áreas A, e Ay isto 6” O Parâmetros são obtidos 


A= | LAy(co) = y(t) dt 


A, = [aytar iaad. 
Hj Ay(o0) 


conforme 
Fist “lustra a figura (4.6), Assim, 
“8 Parâmetros do modelo são d&” 


Ke dy 
Au 


(410) 


(411) 


(412) 
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A, 
Poderoso a 
0.368Ay(00) (4.13) 
O=t;—7 (4.14) 


- Método de Nishikawa para a modelagem de processos de primeira order 


Figura 4,6 


(tb) 


y(0) 


Observações: 

* Por dependerem do traçado da tangente, os métodos de ZNe 
HAG apresentam sensibilidade na presença de ruído (impre- 
cisão na calibração dos parâmetros); 

* Os parâmetros estimados são selecionados de tal forma que 
a soma dos quadrados da diferença entre os valores medidos e 
os calculados é mínima. A função custo (no sentido dos míni- 
mos quadrados) é dada por 


Ide 4.15) 
J=>>.e“(i) (4 
nã 
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na modelagem de um processo de primeira ordem. 


Exemplo 4.1 - De acordo com a curva de reação do processo da figura 
(4.7), estimar um modelo paramétrico de Primeira ordem pelo método 
de Hágglund (admitir uma entrada degrau unitário). 


tes parâmetros: 


SISTEMAS ( 


IDENTIFICAÇÃO DE 


tras do experimento e e(i) r presen, 


onde N é o número de amos bservada (LJUNG; SÓDERSTROM 


o erro entre as saídas real e Ol 
1983); ; sta ao degrau reais, os parâmetros p 
inte Fei poeira pc dependendo das conde 
Ei Papeis do processo, isto é, da magnitude da mu: dança 
do degrau de entrada e da direção da mudança. Estas van 
ções podem ser atribuídas às não linearidades do processo, 


à seguir, o exemplo (4.1) ilustra a aplicação da técnica de Hágelund 


A partir da resposta ao degrau unitário, identificam-se os seguin- 


Ganho 
valor final = 2 3 K=2 
Atraso de transporte 

> 60-12 min 
Constante de tempo 
83.2% do valor final é1.26 


1.26 ocorre em t = 4.35 min 
T=435-1,2 im 
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Jinâmica do processo para a aplicação do método de Hagglund 


- 


O 25 5 75 1 125 15 175 20 225 25 275 30 


tempo (min) 


Exemplo 4.2 - Considere o seguinte processo de quinta ordem repre- 
sentado por: 

277.5 
sº 419.55! +1415º + 4535” +6085+277.5 


G,(s) = 


- O teste de resposta ao degrau unitário fornece uma curva de rea- 
São, conforme ilustra a figura (4.8). ; 
Utilizando-se o primeiro método, de Ziegler/Nichols, pode-se iden- 
tificar um modelo matemático de primeira ordem: K, = 1, 0 = 0.6893 s 
ET=243795. A figura (4.9) compara as respostas dos processos real 
modelo de quinta ordem) e estimada (modelo de primeira ordem) para 


“m sinal de entrada degrau unitário. ; 
Percebe-se que dao identificado não estima Requadame ae 
pi go a 
Jichols são recomendados. Obviamente, há diferentes 
“Normação do processo a partir da curva de reação da resposta ao degrau. 


Digitalizado com CamScanner 


IDENTIFICAÇÃO DE SISTEMAS DINÂMICOS Ly, 
us) 


88 


saidas do sitema e do modelo 


resposta ao degrau unitário 


Figura 4.8 - Resposta ao degrau do sistema 


Figura 4.9 - Comparação das respostas real e aproximada 
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Figura 4.10 - Exemplos de processos de primeira ordem 


válvula controlador 


* processo de nivel 


+V 


amplificador 
de potência 


amplificador 
operacional 


—V engrenagens 
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Alguns exemi 


4.10). 

jos na figura ( , : 

rs aidetando O servomecanismo de velocidade, onde ;, « 
o! 


s x Lá 
btém-se a representação em diagrama de blocos dada na figura (4x) 
o! Ro ) 


plos práticos de sistemas de primeira ordem são am 
s e. 


ou seja, 
siaura 4.11 = Diagrama simplificado do servomecanismo de velocidade 
constante do ganho do tacômetro e 
potenciômetro ampl. potência engrenagens 


o Y, 
e at 


janho do motor DC 
ampl. operacional 


cuja representação matemática é dada por 
Vi(s) a K, 


(4.18) 
8,9) tast1 E 
onde ta 
KKKKR mu 
= aa br: 
a segu 
Ra DR 
. it ' 
Ambos os sistemas apresentados na figura (4-10) são caracê na 
rizados por uma equação diferencial de primeira ordem; ENE 
ça pao de primeira ordem pode sr implementado 5 ETR 
analógica) através amplificado; a 
Mustra a figura (412) ki Pense 3 
função de transferência é calculada por 
R 2 
Eo(s) ” ( R, ) e 
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XE E 
——— — então R;= R»=1 MO, C=10 E; 
Se (105+1) E 


Se pet então Ry=R,=1 MO, C=1 uF; 
(s+1) 


Figura 4.12 - Im 


ementação analógica da equação [4.16] 


A implementação da técnica de estimação de Sundaresan/ 
Krishnaswamy pode ser estabelecida por comandos do softwa- 
re Matlab, (tabela (4.2)), ou seja, considerando a planta de 
quarta ordem expressa por G,(s) = 1/(s + 1)! o modelo FOPDT 
tem a seguinte parametrização: K, = 1, 0=2.15se 7= 2.045. 


Tabela 4.2 - Código em Matlab da estimação da planta de quarta ordem 


mação pelo método de Sundaresan/Krishnaswamy 
tf(L,conv((1 2 1), [1 2 2)))5 
step (gp); 

y (end); 

= interpi (y,t,0.353*kp); 
erpl(y,t,0.853*kp) + 

*t1 - 0.29*t2; 

67*(t2 - tl); 

f(kp, [tau 1), 'inputdelay', teta); 
—9P, gm) ; 


duzem a 

- Os métodos de estimação apresentados ss é suficientemente 

Primeira ordem. Se o modelo de primeira ardem constantes de tempo 
Preciso, um modelo de segunda ordem com duas 
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92 a SS e 
aa stimado. A identif; 
exas conjugadas pode ser és us identificação 


raízes compl iliz: écni io 
ou com es de tempo utilizando técnicas Bráficasé 


de mais do que duas constant 
uma tarefa complexa. 


4.3 MODELAGEM DE PROCESSOS DE SEGUNDA 
ORDEM 


Considere a representação matemática clássica de um sistema de 
segunda ordem subamortecido descrito pela seguinte função de trans. 
ferência (OGATA, 2010): 


o, 


G(9)=—>—— E ral <1i 
P +2la,s+0,? para € (418) 


A resposta do sistema à entrada degrau unitário é 


Dem 
y(D=1- Be “cos(o,Bt—4) (419) 


onde 


Y(t,)=1+e 74 (420) 


e acontece no instante de tempo 


t ==! 5) 
pues: (4.24 
0,B y 
Assim, medindo y 
culado pela equação (4 


(to), o fat ca 
, 20) é à er dé amortecimento, 4; pode se! E, 
determinar a frequênc 


ação (4.21) pode ser utili 


| 


ia natural, «9, 


Bote 
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pru 


re principais especificações de desem 

um sistema de segunda ordem são apresentadas na figura (4 13) e 

em ser diretamente mensuradas Por inspeção da resposta temporal 
a tilizadas na elaboração da correspondente função de transferência, 


ou seja, 


penho da resposta ao degrau 


Figura 4.13 - Resposta ao degrau para um processo de segunda ordem 


A Si =] 


E À 


* tr (tempo de subida); é o intervalo de tempo que o sistema leva 
para ir de zero a 100% do valor final para sistemas sub- 
amortecidos. Para sistemas sobreamortecidos, utiliza-se o in- 
tervalo de 10 a 90%; 


t, (tempo de pico) é o tempo gasto para a saída atingir o pri- 


Meiro pico de overshoot 
(4.22) 
Pio) 1-6" a 
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94 es 
po que o sistema leva para en. 


ão): o tem 
* ts (tempo de acomodação): O ; pes 
trar e permanecer na faixa +2% do valor final (referência) 
4 
Lisias (4.23) 
Co, 


* PO (percentual de overshoot, sobressinal) 


PO(%)= Yet) 00% =100e 7 (4.24) 
y(c0) 


* ess(erro em regime) 
Es =y(00)-y(00)= lim s [Yi(s)-Y(s)] (4.25) 


Portanto, manipulando-se matematicamente algumas destas equa- 
ções, é possível, então, estimar os parâmetros de «» e, consequentemente, 
o modelo de segunda ordem da equação (4.18) (DORE; BISHOP, 1995). 


Observação: 
* O modelo discreto da função de transferência de segunda ordem 
gi=2 0 YO=-aytt-D-ay(t-2)+ 


(sta? É byult-D+bu(t-2) 


tem a forma 


G,(2)= 2'B(e") v(z) 


Ale") U() 
onde A(z!)=1I+az! ras 2, Bql 
entes calculados através de (ROFEI 
VISIOLI, 2009) 


)=by+bz! e os coefid- 
EL et al, 1989, FADAL 


=-2e"%, 
aj=-2e ; a=e ut 


=1-e% 
bo=1-e“-(1+aT) ; b,=e (ek pat) 
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quero 

um segundo procedimento gráfico para a identificação de siste- 
s de segunda ordem descritos por respostas oscilatória ou sobre- 
S recida baseia-se no método apresentado por Mollenkamp (1988) 
curva de reação do processo são identificados três pontos 
intermediários, isto é, tr — tempo para a saída alcançar 15% da mudança 
total final; t> — tempo para a saída alcançar 45% da mudança total final; 
«.- tempo para a saída alcançar 75% da mudança total final. Com base 
nestes instantes de tempo, os parâmetros do modelo paramé- 
de segunda ordem são calculados pelo seguinte algoritmo: 


tt, 
tt 


tnco 


— 0.0805— 5.547(0.475- x)? 
(x- 0.356) 


» 


o EO) =(0.708)(2.811)* se C<1 
£(0)=2.6C-0.60 se C>1 


e) EMO =(0.922)(1.66) 
9 0=t,- AS) 
«o, 


n 


[r2 de. e 
Dis Gty6 1 (somente se a condição 621 é satisfeita) 
j [0 


lo de aplicação do método de 
de segunda ordem. 


M A seguir, apresenta-se um exemp 
ollenkamp na modelagem de um processo 
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e 


E INÂMicas Um 
Ay 


do processo a 
- De acordo com a resposta ! Dave 
a a estimar um modelo paramétrico de segunda pr ER 
oie de Mollenkamp (variação da entrada degrau de 5), Pe 
m 


Figura 4.16 = E xperimentação de um processo para estimação Paramétrica 


so 
“o 
a valor final = 49.5 
47 
“6 
vt) 45 
as 
“3 
“2 
a 
od 25 50 75 100 125 150 175 200 225 250 275 300 
tempo (s) 
15% de 9.5 = 1.42 tie = 508 
45% de 9.5 = 4,28 > tux=80s 
75% de 9.5 = 7.12 tra = 1255 
Vlicando ds intento ta Sis 0 do, nalguns do 
kamp, obtém-se 
K=1,9 
0=272s 
n=527s 
w=19.9s 


e Md = se meetsao mu 
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portanto, a função de transferência de segunda ordem estimada é 
O 

k À K,e 

G,(s) i 


P 


1.9e 27.29 
(mst D(ts+1) (19,95+1)(52.75+1) 


A tabela (4.3) ilustra a programação em Matlab do algoritmo de 
Mollenkamp para O exemplo (4.3). 


| Tabela 4.3 - Código em Matlab do exemplo (4 3] 


ã 


o de um modelo de segunda ordem 
pelo método de Moll 
' valor de t1(15 
valor de t2(45%): '); 

valor de t3(75%): '); 

[(t3-t1); 
=(5.547*(0.475-x)*2))/(x-0.356); 


ra: 


“zeta; 


mo 


2.6*zeta)-0.60; 


T2l(t3-t1)% 
(0.922)*(1.66"zeta); 
t2-(£3/wn); 

(zeta+sgrt (zeta”2-1)) /wn; 
(zeta-sgrt (zeta”2-1))/wni 
Zeta < 1 

disp('fator de amortecimento!) 
zeta 

Cisp('frequência natural!) 

Wn 

disp('atraso de transporte") 
teta 
Else 
Sisp('constante de tempo 1") 
talz 

disp('constante de tempo 2!) 
tal2 


disp('atraso de transporte!) 
teta 
end 
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gráfica alternativa para estimar uma 

em é o método de Smith, onde Sobrea 
curva de reação do processo são identificados dois instantes de te 
t; e t», correspondentes a 20% e 60% do valor final (com o aparente 
atraso de transporte removido). Com o auxílio da figura (4.15), obtém. 
se & com o valor de tzo /teo e a estimativa de 7 é calculada a partir de 
(tso /7) x (too /teo). 


Uma terceira técnica 
de transferência de segunda ord 


Figura 4.15 - Identificação dos parâmetros tzog, toom 7, 6 pelo método de Smith 


Para sistemas sobreamorte, elaDA éguas 
Ke! 
=—— Po. 
(ts +1)(r,s+1) * Ma=t6tr fera 
enquanto que, para sistemas suba; À 
cidos, em-se 


8,9) 


8 (9)= s “o er 
Ts +2ts+1 
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Um exemplo de sistema de controle de segunda ordem de grande 


utilização prática em controle de processos está ilustrado na figura (4.16). 


Figura 4.16 = Sistema de controle de posição analógico clássico 


amplificador 
de potência 


amplificador 


o |] mot e: 
N peracional tor 
engrenagens SE 
+ e VANESSA 
% de posição 
V 
Observações: 


* O sistema servomecanismo de posição está caracterizado em 
malha fechada por uma equação diferencial dominante de se- 
gunda ordem; 

* Aimplementação analógica de um processo de segunda ordem 
através de amplificadores operacionais é apresentada na fi- 
gura (4.17), isto é, 


Figura 4.17 - Exemplo de um processo de segunda ordem 
18MO 
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ati DOS Leg, 


100 o ] q - ; 
j 0.2, 0.3, 1.2) na implementass 
Admitindo op rr ape pdsre qç ão á 
esso de segunda ordem da ; « 
raça do tipo sobreamortecido, criticamente amortecido e sub, 


i amor. 
tecido, Estes resultados estão ilustrados na figura (4.18). 


iqura 4 18 - Respostas lípicas de segunda ordem 


a ganho = 1.2 
saida | 
afetos 
E ro=02 ps 
ganho = 0. presos | 
04 ES do omé 
02 ganho = 0.2 “posta ey 
OIE GEARS 25 25 A 350 am 
tempo (s) Os 
4.4 MODELAGEM DE PROCESSOS VIA RESPOSTA ! 


EM FREQUÊNCIA 


quência é obtida pela substituição de eras mica a resposta 
0) isto é, 
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G,Go)=R(0)+ iX(0) 


(4.28) 
representada em termos da magnitude e fase, ou seja, 
G, (Go) = G, (9) Zó(0) = G, (Jet: ) es 
G, Go) = ÍR(o)P + Ix(o)P is 
[X(0) | 
do)=t, y Fdateod 
e LR(o)] (4.31) 


Em aplicações práticas, deve-se garantir entradas senoidais em 
diferentes frequências e medir precisamente nestas frequências as ra- 
z6es de magnitude e o deslocamento de fase. O método é aplicável para 
processos lineares e o procedimento de identificação é off-line. Além 
disso, o método é aplicável somente para sistemas estáveis, desde que a 
Tesposta em frequência de sistemas instáveis não pode ser medida na 
Prática (DORF; BISHOP, 1995; OGATA, 2010). 


Observação: 
* Um sistema linear e invariante no tempo estável sujeito a uma 
entrada senoidal da forma (u(t) = Usen(at)) apresenta, em re- 
gime permanente, uma saída senoidal, isto é, 


YO = Ysen(ot+ 4); Y=UG, Go) 


do sinal de entrada, porém com a 
se da saída diferentes do sinal de 


(4.32) 


com a mesma frequência 

amplitude e o ângulo de fa: 
entrada. 

ia de 

af de vet a cai da rea em 

“Mas dinâmicos contínuos é necessário ds c/bos UA 


(pol em termos dos traçados de magni 

ze e a função de transferên- 
fia, é € zeros) pra encontrados em uma e 

OU sei 
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tm E Ss Ng, A 

Ganho: K, 

O ganho puro, Ky tem amplitude 20logK, (dB) e fase O. pa à 
qualquer frequência, isto é, transfere igualmente para a saída tod 2 “a 
frequências de entrada e sem defasagem. 

My 
so 
M(dB) 20 E 
Í ar 
| nº 
fr a FE - o 0 
frequência (rad/s) 
=180? 
101 


Polos (ou zeros) na origem: (jay)N 
Admitindo o caso de 


—20log, dB (N20log «) ; 
emw=10 3 -20 
ra 

e a amplitude é uma reta com incliy Gs a 


o) =-90' (90'N) 
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(iu 


10º 10 10º 10º 
frequência (rad/s) 


(jus? 


10 10º 10º 10? 10º 10+ 
frequência (rad/s) 


Polos (ou zeros] sobre o eixo real: A+rjoN 
Considerando o caso de um polo simples (N = -1), então 


—10log(1+w?t?), dB 
HKo)=-tg+or 


V<< 1/7 3S0dB 0=020º 
2>>1/7 3 -20log qr, dB 00500 -90º 
0=1/7 3-3dB 0=1/79-45º 


O módulo proximado por duas 
“vínto do fator polo simples pode ser aí dios 
ncia “ma para baixas frequências de O dB e uma para altas fre- 


de 
“te qa ge inclinação -20 dB/dec com cruzamento na frequência 
T. 
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bad 
frequência (rad/s) a 


WT 10m 


out 


Digitalizado com CamScanner 


-, Métodos clássicos para modelagem de processos 


opiruDs TE es == 105 


ou) 


047 mM 10 100/T 1000 
frequência (rad/s) 


Polos (ou zeros] complexos conjugados: 
5 NaN 

(E os dead 

(to, Da 


Admitindo um par de polos complexos conjugados, N = -1, então 
Grigu-pya, =D 
o, 


20l0g/G(0)]=—10log(1 — 2)? + 40212) 


H<<1 > magnitude = 0 dB, (0) = 0º 

H2>1 > magnitude = - 40 log q, d(t9) = -180º 
BA 

“ou 13 magnitude na vizinhança de O dB 

depende de €;, (mm) = -90º 


P, 
Mage 0.707, a frequência de ressonância é calculada por 


q, =0w,/1-26 
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SISTEMAS DINÂMICOS Lie; 
O Ma 


[Quo 


to que à magnitude máxima de G(jay) é dada por 
enquan 


1 
Max 2i-o 


= 0.8, tem-se a seguinte resposta em frequência: 


Sec 


Ow, IR 1Ousiy, 


frequência (radis) 


frequência (radis) 


Atraso de transporte: e * 


A resposta em frequência do atraso tem módulo unitário * 
fase igual a -«vl), gi 


Exemplo 4.4 - Considere a função de transferência do processo eo o” 
respondente diagrama de Bode traçado na figura (4.19). 


Gs 1.89 
ld)= 5º +5.95! 412.035" +9.405s+1.89 
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3 4.19 - Diagrama de Bode do sistema de quarta ordem 


Wç=1,17 


—80 dB/dec 


frequência (rad/s) 


E E i s frequência de 
Pelo dia “ itude, identificam-se uma 
“orte em 117 radis é Ea de inclinação de -80 dB deco lito comia 
“o seguinte modelo matemático aproximado: 
So (= (0,117 
síntota de alta fre- 
O ganho estático do sistema é e E A inclinação de 


Mência corta emo= 
a reta do ganho estático ão ao número de 
/dec sugere ao sistema quatro a em relaç do modelo esti- 


2eros aproximação 

- Para comparar e ter uma diagramas de magnitude e 
Pr Perante o sistema real, traçam-se 08 

“onjuntamente na figura (4.20). 


| 
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tas em frequência do sistema e do modelo estimado 
ostas « q 


fase (graus) 


frequência (radis) 


Informações adicionais sobre a teoria da resposta em frequência 
podem ser obtidas de forma detalhada nas referências Dorf e Bishop 
(1995), Levine (1996), Faleiros e Yoneyama (2002), Ogata (2010). 


4.5 MODELAGEM DE PROCESSOS VIA RESPOSTA 
IMPULSIVA 


o; métodos de identificação descritos nas seções anteriores pres 

e senoidal, Esta seção apresenta 
parâmetros da função de transferênde 
sta impulsiva (sequência de ponderação) 
c t ser 


Digitalizado com CamScanner 


modelagem de processos 109 


a) As medidas da sequência de ponderação estão mensuradas n 
universo de interesse do experimento; É 


b) O nível de ruído nas informações do processo deve ser míni 
ou seja, sinal/ruído > 5; R 


| c) A ordem do modelo do processo deve ser conhecida a priori. 


| Considere a função de transferência discreta de malha aberta (a ser 
estimada) da planta G , (2) de ordem n e descrita por 

Ye) bo+bz! +b,z? +.tbz” 

U(z) * l+a;z'+az"+..taz” 


(4.33) 


Admitindo uma sequência impulso unitário na entrada do pro- 
cesso (u(t)=1 somente para t=0), então a saída é a resposta impulsiva e 
alculada por 


Be) = Pe! = (0) + he ++ 430 
k=-0 


Com as equações (4.33) e (4.34) é possível derivar a seguinte relação: 
icia bo +b,z)+b,2? +...+ bz” 
á = (0) + [(1) + a, O ++ 


(A Ea a 


[rem Sasn-o fr +... (m)n) 
a A 
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EMAS DINÂMICOS Lys 
o mlNE ge 


Ó0s fre0] h(o)7 
1 O ss (OMI HO) 


(4 38) 


persa: o [he | 


Ei 


Da equação (4.36), observa-se que os coeficientes do Polinômio 
B(z") são identificados se os coeficientes de A(z”) são conhecidos, Para 
determinar os a; (i = 1, ..., n), considere os termos de z"*“ até 2%, Desde 


que não existem outros termos by (i = 1, ..., n) no lado esquerdo ds 
equação (4.35), é possível escrever 
[hem h(2) .. hn) Ta, —h(n+1)] 
hQ) h(3) . h(n+1) |a,, | |-h(n+2) 
e. ... ... E (é E) 
h(n) h(n+1) ... h(2n-1) a a] 


A solução da equação (4.37) existe se o rank da matriz contendo 
sequência de ponderação é igual à ordem do sistema, n. Portanto, 
resolvendo-se primeiramente a equação (4.37), identificam-se os coef- 
cientes do polinômio A(z?) (polos de malha aberta) para depois, com à 
equação (4.36), calcular os coeficientes do polinômio B(z!) (zeros de ma 
lha aberta). 


Exemplo 4.5 — Obter as equações (4.36) e (4.37) para um processo car” 
terizado pela seguinte função de transferência discreta e sequênda e 


ponderação: 
fico) = Potbiz! +b? 
lrajz!ra,z? 
hOz * + h(3)z 3 (aja A 


sê 
De acordo com a equação (4,35), é possível escrever (ignorando 
os termos acima de 2") 


=h(0)+h(1)z ! + 


” 


= 


e 


sésit 


E, 
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quintos "EE agem de processos 


1m 


1 2 F: E 
bo +byz" +b;z? =h(0)+[h(1) +a h(0)b + 
[a,h(0) + ajh(1) +h(2)p ? 


[asha(1) + a, h(2) + 63)? + fa, (2) + a,h(3) + (4) * 
la,h(3) +a, h(4)b +fa (4) 5 


Relacionando os coeficientes de (n + 1 = 3) até (2n = 4), então 


[ee hO)"I-h6)] 
a, h(2) h(3) ] —h(4) | 


Os coeficientes do polinômio B(z”) são obtidos por 


b] [1 0 ofho] 
bl-ja 1 0/hD, 
bo) la a 1/2] 


Exemplo 4.6 - Identificar uma função de transferência discreta e um 
Processo causal de terceira ordem admitindo que as amostras da se- 
quência de ponderação são: h(0) = O, h(1) = 7.157098, NB) À 9.491077, 
h(3) = 8.563889, h(4) = 5.930506, h(5) = 2.845972, MO) =0.144611. 

A função de transferência a ser estimada tem a forma 


ed 2+bz” 
à (0) = Po thus tdos ae 


-3 
Z)= 1 po 
A ) 1+a;z +az +asz 


Da equação (4.37), é possível avaliar os valores de as, as e as, isto é, 


7157039 9.491077 8.563889] * -5.930506 | 
Ago? | 5.930506, 


—2,845972 


cn 


> Ergo 
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Ea ; a,;=-04965 
a/==2082575/3) Ho =.764089/ às -0.496585 


Substituindo-se estes valores na equação (4.36), obtém-se q A, 


guinte polinômio para B(z 1 
B(z1)=7.1570312” — 6.487547z * 
Neste caso, a função de transferência discreta estimada é dada por 
a 71570312" - 6,4875472 * 
(= 22325752! +1.7640882? 0,4965857” 


A tabela (4.4) ilustra a programação em Matlab do algoritmo de 
estimação da resposta impulsiva para o exemplo (4.6). 


Tabela 4.4 - Código em Matlab do exemplo [4.6] 


* Estimação pelo método da resposta impulsiva 
$ do exemplo (4.6) 
na=3; nb=4; 
h(1)=7.157039;h(2)=9.491077;h(3)=8.563889; 
h(4)=5.930506;h (5)=2.845972;h(6)=0.144611; 
for i=l:na 
for j=l:na 
ha(i,5)=h(i+j-1); 
end 
hauxl (1)=-h(i+3); 
end 
=i f =] * te 
rt cb ahchal*hauxl'; a=[ah(3) ah(2) ah(1)] 
for i=1:nb 
if j>i 
hb(j,i)=a(j-i); 
elseif j==i 
hb(j,1)=1; 
else 
hb(9,1)=0; 
end 
end 
end p 
haux2=[0 h(1) h42) h(3)7; 
b=(hb*haux2')! 
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|) 
Na cupiruLDA M 
46 PRO BLEMAS 
diferentes métodos de estimação de modelos paramétricos 
rdem na identificação dos seguintes processos: 
A 
e 


SI * Gy(s) 


| naplicaros 
de primeira o 


“(+ 


2) Considere o sistema de terceira ordem representado por 
10 
s(0.2s + 1)(0.5s +1) 


jonar em termos de módulo e fase e obter os valores exatos de am- 
plitude e fase para a faixa de frequência de 0.5 a 50 rad/s. Plotar as curvas de 


resposta em frequência e comparar os resultados exato e assintótico. 


3) Obter a função de transferência estimada para as seguintes curvas de 


Gols) = 


magnitude: 
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114 si 
co) 


“0 


M(dB) 20 


10* 
frequência (rad/s) 


4) Considere as medidas da resposta em frequência para uma combi. | 
nação de servomotor e servoamplificador em um laboratório de controle | 
de processos, conforme apresentado a seguir. | 


Plotar as curvas da resposta em frequência a partir dos dados da tabela 
e estimar uma função de transferência para o processo. 
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30 -3.7 
4.0 -6.0 lis 
ss 5.0 =7,9 
| 
7.0 =10.8 
10.0 | =138 1) 
Coy 5) As amostras da resposta em frequência são apresentadas na tabela. 
atm Considerando a seguinte função de transferência da planta: 
K.(s+z 
e p(8+2) 


“(+ +Zo,s+o?) 


plotar as curvas de Bode e identificar os parâmetros K,, z, p, & «» do 
sistema. 
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) Quando uma entrada degrau foi aplicada em um sistema de segun. 
6) Quai 


da ordem, a resposta em regime foi medida em 2.85. A sobre-elevação 
a O , à TeSpoS 
máxima foi de 20% e ocorreu em 25 ms após a aplicação do degrau. Iden. 


tificar o modelo do sistema. 


7) As amostras da resposta impulsiva de um processo de segunda ordem são 


Rd o joos[01 ne 0.2 a 0.3 [ass 04 
MO] o [068[107[1.27/134/133/128[121 [111 


Determinar as funções de transferência contínua e discreta da Planta, 


8) Um engenheiro de controle de processos de uma fábrica de motores 
realiza um ensaio pela resposta ao degrau de malha aberta em um ser- 
vomecanismo de posição. Aplicando um degrau na entrada do sistema 
(tensão de alimentação), observa a resposta conforme a figura (4.21). 


Figura 4,21 - Resposta ao degrau unitário de malha aberta 


vít) 
RE Rs 
ut, y(t) 
gs [ sistema ]25 t 


Figura 4,22 - Resposta ao degrau unitário de malha fechada 


Não identificando os parâmetros do sistema, o engenheiro decide 1” 
Plementar um sistema de controle proporcional, conforme luste * 
figura (4.22), realizando sucessivos ensaios pela resposta ao degra! 
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E" fechada com diferentes valores de K.. A, 
m o avaliar o modelo matemático do sist 


Pós várias respostas, de- 


8 qíde entã! ema como ilustra a figura 


k (422) onde K. = 0.5. 
a) Qual o modelo matemático estimado obtido pelo engenheiro de con- 
trole de processos?; 
b) Mostrar que O ganho estático do sistema é unitário e o erro em regi- 
me é nulo; 
carga do sistema de posição de malha aberta para uma 8 
entrada pulso unitário com duração de 10 s. 


| 9) Um engenheiro de controle realiza diferentes ensaios sobre um sis- 
tema de controle. Na primeira experimentação, aplica um degrau uni- 
tário. Como a resposta é ruidosa, o engenheiro somente consegue deter- 
minar o valor final da resposta. Neste caso y(%) = 2 e a derivada na origem 
de y(t) é nula. Na segunda e na terceira experimentação, obtém os se- 
Buintes valores em regime: 


ER 


u,(t)=sent — Iya(o) = = 


N 


A 


u,(t)=sen?t > |y,(t)] E 


J10 


minar a função de transferência contínua do sistema a partir 
dp informações dos três ensaios; 


e . 
Pi 9 sistema, implementa-se uma lei de controle PI digi- 


3 


—6 YB=u(t-1)+go[e(t)=0.905e(t=)] 
& | Mismii MP=IA-1O 
ego 


ses | E 2 € (0) é a sequência de referência do sistema. Ajustar O 
vs! | Praqueo sitema de malha fechada apresente estabilidade 


7” 
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a FR RE - TO aa 
engei cessos obtém sete 
nheiro de controle de pro eis 
10) Um vimental com um filtro digital exponencia] E 4 


meira ordem em uma planta de temperatura de acordo com os da dt 


seguinte tabela: 


A 


5 JR pEM Sn 
1644 | 1575 | 1533 | 1531 ERR 
ú 


167.3 | 165.3 | 1629 | 160.9 | 1601 [igoy 


Considerando y(t) a medida real, y(t) a medida filtrada e T, = 05 mao 
período de amostragem, então: 

a) Qual é a frequência de corte do filtro aplicada na experimentação?. 
b) Qual é a medida que recebe a maior ponderação, y(t) ou y(t)? 


11) Seja o filtro passivo RCR 


Ro 
Tê 


vo] [o 
R2 


onde u(t) é a entrada e y(t) é a 
É saída. Co; em 
a) Identifi obtida em ensaio de laboratgmoo e Pá SeBuinte peer 
aestrutura ea a 
Parametrização da função de transferênd 
filtro, incluindo os valores de R7 e R$ considerando po 
b) Qual é a característica do filtro? ? 


O mn cn A 
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pot 


107 10º 10º 10º 10º 
frequência (rad/s) 
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CAPÍTULO 5 


IDENTIFICAÇÃO 
DE SISTEMAS 
REPRESENTADOS 
POR EQUAÇÕES 
A DIFERENÇAS 


ã delos 

- Diferentes métodos timação de parâmetros de mo: 
para a es | 

qts discretos são apresentados. Ênfase maior é dada ao estimador 
de, Mimos quadrados, uma vez que é a base para o desenvolvimento 
“tros métodos de identificação (LJUNG, 1999). rins 
ey A Qualidade da estimação é dependente da natureza á - sr 
“tura do modelo, do tipo de aplicação e da “riqueza” da informa: 
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contida nas medidas. Tanto a implementação off-line quanto à on-line dy 
ectos computacionais são avaliados 


UEN, 1990; BOBÁL et al. 2005), 
vam a utilização do estimador 
) no teste da resposta ” 


estimador são apresentadas eos asp 
(ISERMANN, 1980; HABER; UNBEHA! 

Algumas características que moti 
dos mínimos quadrados são as seguintes: à À 
degrau, as medidas do processo são obtidas a partir da mudança de um 
ponto de operação para um novo ponto de operação. Como em apli- 
cações práticas têm-se processos não lineares, O deslocamento do ponto 
de operação desejado é perigoso do ponto de vista da estabilidade, e 
também pode proporcionar produtos manufaturados fora das especi- 
ficações de projeto; b) o método do teste da resposta ao degrau é ade. 
quado se o processo sob avaliação contém um baixo nível de ruído 
(variância reduzida), o que muitas vezes não acontece em condições prá- 
ticas industriais; c) de forma diferente do teste de resposta ao degrau, 
cuja interpretação gráfica é aplicada no conjunto de amostras coletadas 
da experimentação (procedimento para inferir sobre os parâmetros do 
modelo, onde não necessariamente a magnitude da variável do processo 
medida “domina” o nível de ruído), o estimador dos mínimos quadrados 
manipula as medidas de entrada e saída nas formas não iterativa e 
iterativa através de algoritmos não recursivo e recursivo. 

Considera-se a priori que a ordem do modelo é conhecida e que as 
amostras das medidas de entrada e saída estão disponíveis a cada perio- 
do de amostragem no universo da experimentação (equispaced). 

Os principais objetivos deste capítulo são: a) familiarizar o leitor 
com a derivação, propriedades e utilização do algoritmo dos mínimos 
quadrados recursivo (MQR) quando aplicado na estimação de sistemas 
lineares discretos SISO (Single-Input and Single-Output) com parâmetros 
desconhecidos e/ou variantes no tempo; b) na impossibilidade de j 
cação do estimador dos MQR, estudar algoritmos alternativos de est 
mação paramétricos; c) viabilizar a utilização do estimador dos MQR e” 
aplicações de controle adaptativo, isto é, no contexto de algoritmos 
controle autoajustável (selftuning control) (ROFEBL et al., 1989; ASTRÔM: 
WITTENMARK, 1995; COELHO et al, 1999), 
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«2 FORMALISMO HISTÓRICO DOS MÍNIMOS 
” QUADRADOS 


Karl Friedrich Gauss formulou o Princípio dos Mínimos Quadrados 
ao final do século 18 para prever a trajetória de Planetas e cometas a 
partir das observações realizadas. K. F. Gauss estabeleceu que os pará- 
metros desconhecidos de um modelo matemático deveriam ser selecio- | 
pados de modo que a 


o valor mais provável das grandezas desconhecidas é a que mini- 
miza a soma dos quadrados da diferença entre os valores atual- 
mente observados e os valores calculados multiplicados por nu- 
meros que medem o grau de precisão, onde quanto mais precisa a 
medida, maior a sua ponderação. (LJUNG; SÓDERSTRÔM, 1983). 


5.3 ESTIMADOR DOS MÍNIMOS QUADRADOS 
NÃO RECURSIVO 


Considere um processo físico caracterizado por uma entrada, u(t), 
uma saída, y(t), uma perturbação, e(t), e com função de transferência 
discreta linear na forma 


A(z y(t)=2 "Blz ul+e(t) (5.1) 


onde 
a 


A(zl)=1+a;z! +... ta 


Blz )=b,+b;27+...+ [ps 


“cuja representação por uma equação a diferenças é 


y(b=-a,y(t-1)-a,y(t—2)=..— asay(t—na)+ (5.2) 
bu(t-d) + bu(t-d=1)+...+ bpult=d = nb) +e() 
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Observações: 

* Para o modelo da e 
metros a estimar; 

» Para determinar os ai (i= 1, «=» na) e by (j=0, ... nb), deve-se 
utilizar as medidas de entrada e saída do processo; 

* O termo e(t) pode representar o erro de modelagem, o erro de 
medição ou o ruído na saída do tipo estocástico, determi 
nístico ou offset. 


quação (5.2), têm-se (p = na + nb + 1) parg 


Definindo-se o vetor de medidas, «Xt), de dimensão (na + nb + 1)x1 
0 (D=Ly(t-D)-y(t-2)...-y(t-na) ult-d)...u(t-d=nb)] (53) 
e o vetor de parâmetros, G(t), de dimensão (na + nb + 1)x1 
q'(D=[a,a,...a,, Do D,...b,y] 
pode-se reescrever a equação (5.2) como 
y(D= q (DO +e(t) 


que K denominado modelo de regressão linear (LJUNG; SÓDERSTRÔM, 
1983). 


Admitindo que são realizadas N medidas, suficie: u er 
minar os parâmetros a; e by, então tem-se que ; ai 


v(o) 9"(0) e(o) 
y(1) | 9) vá e(1) 
YN-D] |p"(N-1) “N-1)] 
A 
representação matricial da equação (5.6) é 


Y=40+B 
onde a matriz de observação é 
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[-ren  -yta = y(-na) u(-d) ul=d=1) u(-d=nb) 
| yo) 61) vma)  u(l-d) u(-d) u(1-d=nb) 
+=| xD v(o) vma) ud) u(i=d) o u(2-donb) 
| y(N=2 -y(N-3) =Y(N na -1) u(N-d-1) uN-d=2) u(N-d-nb=1) 


e o vetor de saída é dado por 
Yº =[y(0) y(1) (2)... y(N=1)] 
Observação: 


* A matriz $ não é quadrada, isto é, tem mais linhas do que 
colunas. 


A estimativa do vetor de parâmetros, À, pode ser obtida pelo pro- 
cedimento dos mínimos quadrados (least squares approach). Utilizan- 
do a estimativa (), a melhor predição da saída do sistema, 7, é calculada por 


Y=40 (5.8) 
eo erro de predição, &, é avaliado de acordo com 


e=Y-Y=Y-40 (5.9) 


O estimador dos mínimos quadrados ponderado, também deno- 
minado estimador de Markov, é obtido minimizando o seguinte critério: 


empe-4i ão 


J=[7 40] WIY 40] 


ou 


(511) 


e sendo a matriz W simétrica e definida positiva, isto é, 
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fw(o) O 0 
e o w(l) 0 
o 0 wN-D] 18 
Ed 
onente do erro e função da py 


de w() é a ponderação em cada componei 
lodo da medida (quanto mais precisa à medida, maior é a ponderação). 


Derivando a equação (5.11) em relação aú e igualando-a a zero, 
tem-se que 
=a(vTwo) +26" W$0=0 


Assim, o estimador dos mínimos quadrados ponderado é 
por 
O=[4 Wo] 4TWy 
e isso conduz ao mínimo se 
os? 
E = 24" Wy >0 
e, portanto, (4"W%) deve ser uma matriz não singular (full rank) 
Observações: 
* I=W-9) We 0d) = 
T A À a 
a WE VW Dr WY + DrTwad, 
- W=w', 


Considere os vetores b 
EL o 
) mostrar que (M 
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custru 
dor 
——(m Am) =2Am 
dm 
doar 
—(b'm)=b; 
Sr m)=b; 


» O vetor Ô pode ser obtido desde que (4'W) seja uma ma- 
triz definida positiva. A inexistência desta condição está 
associada à presença de medidas insuficientes, incorretas, 
redundantes ou à falta de uma excitação adequada no sistema 
(LJUNG, 1999); 


* O sistema deve estar persistentemente excitado para evitar O 
caso de linhas comuns na matriz é (colunas linearmente de- 
pendentes). 


O estimador dos mínimos quadrados não recursivo é obtido admi- 
tindo que 
W=Lan 


isto é, a mesma ponderação é aplicada em todos os erros de medida. 
Logo, a equação (5.12) torna-se 


ô-[5 fra elo 
O=[4"0]"4"Y (5.13) 


O estimador dos mínimos quadrados, equação (5.13), é uma trans- 
formação linear sobre Y (função linear das medidas) e, assim, é deno- 
minado estimador linear (GOODWIN; PAYNE, 1977; LJUNG, 1999). A 
tabela (5.1) ilustra um código, em Matlab, da implementação do esti- 
mador dos mínimos quadrados não recursivo aplicado na modelagem de 
uma planta discreta de segunda ordem. 
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Tabela 5.1 Programação em Matlab do estimador dos mínimos quadrados: 
ela x 
3 Esciasdor dos mínimos quadrados não recursivo 
dem & 

ocesso de segunda Or y ' 
is medidas contidas no arquivo medidas .dat 
clear all 
joad medidas .dat 
npts=30; u=medidas (1:nptS, 1);y=medidas (L:npts,2); 
ye[J;fi=[); 
for j=l:npts 

1E j<=2 

y1=0;y2=0; ul=0;u2=0; 
else vl=y(3-1) :y2=y (9-2) ul=u(j-1) puZ=u (5-2); 

end; 

Ea y(9));fi=[fi;=y1 -y2 ul u2); 
eng; 
teta=inv(fi!*€i)*fi'xy 
for t=1:2, 

yest(t)=0; 
end; 


al=teta(1);a2=teta (2) ;bl=: 
FR RS a(2);bl=teta(3);b2=teta (4); 


vest (t)=-al* Es 
Sd altyest (t-1) -aZtyest (t-2) 4bl tu (t-1) +b2*u(t-2); 


Plot (y, 'g'); 
hold on; 
Plot (yest, 'r'); 


Observações: 


* A condição de 
condição de excitação," “Niz (FA) é inversível é deno 
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Figura 5.1 - Procedimento geral para a identificação de processos 
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A dimensão de (F') é 
pxp= (na +nb+1)x(na+nb+1); 


Uma vez calculada (estimada) a função de transferência dig. 
creta, a determinação da correspondente função de t 
rência contínua pode ser obtida por uma das seguintes relações 


z=I+Tys — aproximação retangular backward 
z=1M1-Tys) — aproximação retangular forward 
z=(1+T,s/2)M1-T;s/2) > aproximação trapezoidal 


As principais etapas envolvidas na identificação de sist 
dinâmicos estão ilustradas na figura (5.1). Os passos 
nados no diagrama de identificação são abordados neste 
tulo, juntamente com outras técnicas aplicadas na estin 
de parâmetros (LJUNG; GLAD, 1994; BROSILOW: JOSE 
2002; MORAES, 2004); 


Na aplicação do estimador dos mínimos quadrados, figura 
todas as medidas devem estar disponíveis a priori para an 


e não existe limitação no tempo de pr: amen! 
tmo (identificação off-line). dita es ea 


Figura 5.2 - Procedimento de identificação off-line 


R 
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A seguir, os exemplos (5.1) e (5.2) discutem algumas caracterís- 
ticas importantes do estimador dos mínimos quadrados dado pela equa- 


ção (5.13). 
Exemplo 5.1 - Seja um processo sem atraso de transporte representado 
pelo modelo FIR da forma 


y(O=bou(D+b,u(t-1)+e(t) 
A matriz $ para N amostras é caracterizada por 


u(t) u(t—1) 
: Ey u(t+1) u(t) 


u(t+N) u(t+N- 1)) 


para os parâmetros estimados dy e by. 


A matriz (fd) é então calculada por 
u(t) u(t-1) 
“4 | u(t) u(t+1)  u(t+N) | u(t+1) uít) 
u(t-1) u(t) .. u(t+N-2)] 


u(t+N) u(t+N— 9] 


tN Apa E o: ] 
” Zu (1) Qu stDu(d=1) 


i 
tan 


“En 
Duli-Iuli) Su”(i-1) 
+= tt 


Observações: 

* A dimensão de (4) depende do número de parâmetros 
desconhecidos (não do número de amostras). Neste caso, a 
matriz é (2x2). Para m parâmetros desconhecidos, a matriz 


tem dimensão (mxm); 
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:a de modo que somente a pa 


imétrica, : 
Ê (ou inferior) precisa ser calculada; 


ut) = to, então 


“A matriz (P é 
triangular superio! 
Se u(t) é uma constante, 


y(0=(Do +b)uo +e(t) 


a matriz (9) é dada por 
Rea 
po=Nody 1] 


e, assim, a matriz (99) é singular e uma solução úni 
(mínimos quadrados) não pode ser obtida porque (f9) 
inversível. Para a solução ser obtida, então u(t) deve v 
suficientemente para garantir que 


det (676) %0 


quando N cresce, Esta condição está usualmente associada com 
ame peiden ficientemente excitado. No contexto de controle 

, É importante que u(t) mude suficientemente à 
um rank deficiente para a matriz (fd); "4 
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Ma o. E aaa ses 


a 
um simples teste com base no de, 
cado (RIVERA; FLORES, 2000); 


* A precisão das estimativas está associ 
elementos da matriz de covariância, qui 
lada por 


grau ou pulso pode ser apli- 


ada à magnitude dos 
e por definição é calcu- 


PCN) =[4"(N) (ND)? 


* Variabilidade dos parâmetros estimados + 


Exemplo 5.2 - Obter a função de transferência estimada para um pro- 
cesso FIR representado pelo seguinte modelo matemático: 


y(D=bou(t+b,u(t-1) 


utilizando o estimador dos mínimos quadrados (procedimento off-line). 
Determinar as estimativas para Ô dos parâmetros 5, e b, paraN=7,8,. 

14 considerando as medidas de entrada e saída de doa com a do 
ou seja, 
Ejoji[2/sJalsJel7ls 
a o es os tEsfes teste 02 


[0.9/1.4/1.9[2.3[24[2.3/13/1.2 


Para N = 7, tem-se 


[y(D] [u(t) u(0)] 
y(2) | ul?) u(1), 
ya, u(3) u(?) bo 
y(4) = u(4) u(3) [| 
y(5) | | u(s) u(4) 
y(6) | u(6) u(s) | 
[y(7] u(7) (6) | 


Substituindo-se os valores numéricos da tabela, obtém-se 


[a 
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TUR RR 
if popios 1] 
24 0.6 0.8, 
13, 04 do [Bs 
LH 0.2 0. Ê 
08 |0 02 
o | [02 o | 
09] [04 02] 
A estimativa ótima para by e b, é avaliada por 
d=[9"9]20"Y 
e calculada de acordo com Ná é 
0.6 0.8 
Roda mo do O RE qEa DES 
“D=|1 08 08 0. RE Co 04 
(40) Es 08 06 04 02 0 02 EE 
02 0 
04 02 
EE 
cola 168] ; (691 = 7.143 -5357] , 
a -5357 4464 | 
[25] 
24. 
=| 7143 -s357 a 
Lim same | 08 04 02 0 02 04,2 
1 08 06 04 02 0 02o8 
o 
[o3) 
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[0.322] 
2.446 | 


A tabela (5.2) mostra os valores estimados de bp é b, quando se 
aumenta o número de medidas (N). Uma solução por computador via 
Matlab pode ser implementada para parametrizar este problema. Adi- 
cionalmente, a precisão das estimativas de b, e b, pode ser avaliada cal- 
culando-se o traço de P(t) (t,[P(t)]) para cada valor de N da tabela (5.2). 


Tabela 5.2 - Estimativas de dy e b, 


CE E 
7 0.322 2.446 11.607 
8 0.607 2.256 | 7.738 
9 0.661 2.222 5.972 
10 0.573 2.272 5.152 
11 0.662 | 2.116 4.752 
0.606, 2.199 4.483 

2.104 4.246 


2.181 4.001 


5.3.1 PROPRIEDADES DO ESTIMADOR DOS 
MÍNIMOS QUADRADOS 


O estimador dos mínimos quadrados Ô é uma variável aleatória, 


onde as propriedades podem ser analisadas utilizando-se a equação a 
diferenças do processo, ou seja, 


y(D= q" (DO) +e(t) 


a qual define a saída real e a perturbação (ASTRÔM, 1970; LJUNG; 
SÓDERSTRÔM, 1983). As propriedades mais importantes são polariza- 
são (bias) e covariância (covariance), isto é, 
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* O estimador é não - . 
vergem para 08 parâmetros verdadeiros quando o número d 


iterações aumenta) se à perturbação é po ruído io média 
nula e variância 02 e yít), u(t) são estatisticamente indepen.. 


dentes de e(t); 


* Cov(0)= 020)” Moo! “48 
fornece a medida direta da variabilidade e covariabilidade dos 
a precisão das estimativas é estabele- 


parâmetros estimados ( À 
cida pelo tamanho dos elementos da matriz de covariância); 


* Seo estimador dos mínimos quadrados é não polarizado e con- 
sistente (BLUE - Best Linear Unbiased Estimate), então 
zI(ô-00-07] < =I(9*-0)0*-0)] 


(melhor estimativa) (qualquer estimativa) 


E 


=|po-s|sfe-n-e)<...</ôç)-0 


st21 


onde =.) denota o operado! F 
1970). operador esperança matemática (ASTRÓM, 


5.4 ESTIMADOR DOS MÍNIM 
ADO OS QUADRADOS 


A aplicação de 


” algoritm: Cs Ped 
vários propósitos, os de identificação é interessante pata 


tais i 
tomo supervisão, rastreamento de p 


mulação adequada dos al 
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novas medidas estão disponíveis (adequar-se às características do pro- 
cesso controlado ou ressintonizar-se caso existam variações na dinâmica 
do processo). 

Em particular, é importante visualizar o procedimento de estima- 
ção recursivo em termos de um modelo paralelo conforme mostra a fi- 
gura (5.3). 


Figura 5.3 "rocedimento interativo na estimação de parâmetros 
processo 
t 
8 y(t) 
+v 
u(t) 
e(t) 


erro de 
modelagem 


VD = (Mt) 


correção das 
estimativas 


mecanismo de 
adaptação 


Em cada período de amostragem, novas medidas tornam-se dis- 
poníveis e são utilizadas com o modelo atual para gerar um novo erro 
de modelagem (new fitting error), «(t) (mede a qualidade do modelo 
estimado). Por exemplo, no instante de tempo (t + 1), novas medidas 
u(t + 1), y(t + 1) ocorrem. Em vez de recalcular o estimador dos mínimos 
quadrados, é interessante atualizar as estimativas anteriormente calcu- 
ladas no instante t, ó(1), para obter as novas estimativas 0(1=1). 

Observações: 

* A quantidade de dados armazenada no estimador recursivo é 

pequena se comparada com o estimador não recursivo; 


* Outros nomes atribuídos na literatura ao estimador recursivo 
são: estimação sequencial, identificação em tempo real, iden- 
tificação on-line. 
Benim 
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di hor ati iva de uma con: 
ostra à estimativa recurs nm 
O exemplo (5.3) Mm 
i bservação ruidosa de um 
idere a seguinte O à 
Exemplo 5.3 - Consh 


tro constante: y() = 0+ e(t) 


onde e(t) é um ruído branco com média nula e variância o/, À equação 


de regressão linear é 
y(t) = (0 + e(t) 


sendo (e) = 1 (7't >1). A partir da equação (5.13) do estimador dos Mg, 
obtém-se 9/1) como a média das amostras 


a 1Ã..,. ! = 
t)= É » y(i) (614 ti+D= 
1-1 
Para evitar o somatório da saída no cálculo de Ô a cada instantede 
tempo, quando novos dados estão disponíveis, é possível reescrever à 
equação (5.14) como 
Aestim 


d0=0e-D+yt- de) 


a agora 0(1) é atualizado em cada período de amostragem a 
s da saída, y(t) (estimação recursiva). 


3.41 DERIVAÇÃO DO ESTIMADOR MaGR 


Pis º desenvolvimento 


Ê ar a 
variam de 1 (um) Pp ma baseada nas medidas em instantes 


“nstantes de 1 (um) até (t+ Do à estimativa baseada nas A 
mostrado 
ados para amostras pa o (5,3), o estimador dos mínimos y 
deN= =>] é calculado por ( 


de) =[4" E Ty 


Digitalizado com CamScanner 


cpiruLOS Identificação de sistemas representados por equações a diferenças 139 


onde 
9"(1)] y(1) 
dO= q'(2) Ea y(2) 
q"(t) ví) 


Supor que no instante (t + 1) obtém-se nova medida do sistema, 
então os vetores de medida e saída são reescritos como 


q") (1) 
9"(2) (2) e 
t 
dt+D= Ro dE PYQ+D=| .. | et 
s 
O) vo | 4 
q! (t+1) y(t+1) 


As estimativas no instante de tempo t são 


de =[4" (DMD DO 
enquanto que, no instante (t+ 1), são dadas por 
dt + 1) = [47 (t+ DEAD] 4 + DYCE+1) (5.15) 
onde 
dt ] 
q" (t+1)] 
4a D+) = 6" (DMD + o(t+ De" (t+) (516) 


Uma vez conhecido q(t+1), pode-se atualizar a matriz anterior 
das correlações «(t)b(t) para obter a matriz atual (t+1)d(t+1). En- 

to, é necessário encontrar uma maneira de atualizar a inversa de 
F(t) (t) sem calcular a matriz inversa em cada instante de tempo. 


des DAt+D=[6(0 g(t+ v] 
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Adicionalmente, necessita-se atualiza! 


| Y(t) 
grave O MM ye) 
ye (t+ Dy(t+ 1) (54 


ro termo 6(t+1)Y(t+1),istog, 


6 (ta D+) = qr (ore 
Sejam as seguintes definições: 
P(O=[6" (DON 


Rey = 4"(0Y(D 
Substituindo as equações (5.18) e (5.19) na equação (5.15), obtém-se 
d+ 1)=P(t+I)R(E+ TD) 


ou então 


B(t) =PCOR(S) 
enquanto que, das equações (5.16) e (5.17), encontram-se 


PALA D)=PHO+p(t+ Dq" (t+1) 


Rt+D=R(b +q(t+1y(t +1) 

A equação (5.22) fornece uma : 

Dera pi ideia e ren ) 
cando-se a seguinte identidade: pela equação (5.21) 

(ABCD)! =A? - A“B(CA + DASBJ IDA? 

Considerando o termo 
a do lado direito da equação (5:21), isto é 
[POr q(t+1) PACOTE 


onde À = PY(t), C = 
»C=1,B=q(t+1) De aca id) ontémde 


P(t+1)=P(B 
(OU q(t+ (14 qi(t+ DP(Op(t+ 1) grs EO) 
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pÍTULOS. Identificação de sistemas representados por equações à diferenças 
CAPÍTULOS | 
pólo 


tão 
ps PORCA Do" (+ DP 


1+q' (t+ DPCOp(t+1) 


P(t+1)=P(t)— (5.23) 


onde (1 + (t+ DP(D (t+ 1)) é um escalar. 
De acordo com a equação do erro de predição 
e(t+D=y(t+1)-q" (t+ Do(t) 
a equação (5.22) torna-se 
Rt+1)=R(D+ (t+ 1)fe(t+ 1) +qT+ DO) 
R(t+ D=R(D + q(t+ Delt+ 1) + (t+ 1) q” (t+ 1) 6) 


Substituindo as equações (5.16) e (5.20) na equação (5.24), resulta 


(5.24) 


PHt+ D+ 1)=PH0ÔO) + q(t+ 1) (t+ 1)+ 
(P(t+1)-P2(0)Ô() 
Ot +1)= Ô(t) + P(t+1) (t+) e(t+ 1) 


O termo P(t+1)q(t+1) é um vetor coluna e é denominado ganho do 
estimador, ou seja, 


P(tp(t+1) 
1+q"(t+DP(Do(t+1) 


O vetor de parâmetros estimados é calculado por 


K(t+1)=P(t+1)q(t+1) = 


des 1)= O) +K(t+ De(t+ 1) 


nÊ “A equação que calcula P(t + 1) fornece a atualização de P(t) para P(t+1) 
“em inversão de matriz. A única inversão está associada ao escalar 


(+ q" (t+ 1)P(Op(t+ 1) 


4 
p ar Bam 
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A equação recursiva para Pp(t + 1) pode ser combinada com a equ 
ção para Re + 1) de diferentes maneiras para estabelecer uma forr 


i a). ” 
in utiliza-se a definição da variável erro dt + 1) 


como 


eta D=y(t+ Dq" (+DO 


Combinando as equações P(t + 1), R(t + 1) e e(t + 1), obtém-se a 
forma recursiva de Ot + 1) dada por 


Mes) = 0 +K(+D et +D-g' (t+ DÃO) 


á P(Dq(t+1) 
DA DPOt+) 


oro POPt+DOT (t+ DC) 
Met D=O = Cs DPOC) 


A seguir, descreve-se o algoritmo básico do estimador dos À 
a) Medir a saída e entrada do sistema; 
b) Atualizar o vetor de medidas 


(t+ D)=[-y(D-y(t-1)...ult-d+Du(t-d)..]; 
é) Calcular o erro de predição 
t+ D=y(t+ Dq" (t+); 
d) Calcular o ganho do estimador 


Kt+1)= Pole +1) 7, 
1+ q (+ DPOC + 1)" 
e) Calcular o vetor de parâmetros estimados 
de+D=00 +KG A De(e + 1); 
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capirutos Identificaçã 
f) Calcular a matriz de covariância 
POC + Do (t+ DP(O 
ou então utilizar a equação 
P(t+1)=P()-K(t+ IPO (t+ DT. 


P(t+1) =P(t)— 


A tabela (5.3) ilustra a programação em Matlab do estimador dos 


mínimos quadrados recursivo. 
Tabela 5.3 - Código em Matlab do estimador dos MOR 


o de sistemas representados por equações a diferenças 
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YEstimador dos minimos quadrados recursivo 
$ aplicado a uma planta de segunda ordem 
clear all 
nit=input (" Quantas iterações ? '); 
for i=l:nit % Entrada e ruído 
if rand>0.5 
uti)=1; 
else 
u(ij="1; 
end 
end 
e=u*0.01; 
p=1000*eye(4,4);teta=[0;0;0;0]; % Condições iniciais 
for t=1:4 
y(t)=0;erro(t)=0; 


al(t)=teta (1) ;a2 (t)=teta (2) ;b0 (t)=teta (3) ;bl (t)=teta (4); 


end 


for t=4:nit % Iteração da simulação 


Y(t)=1.5144*y (t-1)-0.5506*y (t-2) +0.599*u(t=1)+... 
0.163*u(t-2)+e(t); 
Fi=[-y(t-1);-y (t-2);u(t-1);u(t-2)]; 


erro (t)=y(t)-teta'*fi; % Erro de estimação 
k=p*fi/(1+fi'*prfi); % Ganho do estimador 
teta=teta+k*erro(t); $ Vetor de parâmetros 
P=(p-k*fi '*p) ; 3 Matriz de covariância 
E al(t)=teta (1) ;a2 (t)=teta (2) ;b0 (t)=teta (3) ;bl(t)=teta (4); 
telinit; % Parâmetros estimados 


“bplot (221) plot (t,a? (t)), ylabel ('al') ,xlabel (“amostras”); 
“bplot (222) ,plot (t,a2(t)), ylabel ('a2!) xlabel ('amostras") 7 
lot (223) plot (t,b0(t)) ,ylabel ('b0') xlabel ('amostras 


SMoplot (224) plot (t,bi (t)) , ylabel ('b1') ,xlabel ('amostras"); 


ir 
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tia 
Demonstrar O MATRIX INVERSION L, EMMA 
+ DA!B] matrizes quadradas não sin 


a seguinte relação é válida: 


Exemplo 5.4 


Considere A, G, [C” 
(de dimensões quinas então 


- AÍBIC! +DA 'BJ' DA” 


Bllares 


[A +BCD]" = 

Prova: Pela substituição direta, tem-se 
[A + BCDJ[A? = A “BIC 1, DAÍBY'DA”]= 
I+BCDA? -BCCH(C” + DA BJ" DA? -BCDAB(C" +DABY'DA”. 
[+BCDA? -BC(C! +DA“BHC” + DA“BJ'DA! = 
I+BCDA” -BCDA” =I 

Exemplo 5.5 - Demonstrar a correlação entre o ganho do estimador e 
matriz de covariância. 


K(t+1) = P(t+ Do(t+1) 


E POOq(t+ Do! (t+ DP(t) 
K(t+1)= (poo-Homa + nei) 
neo + nro [E 


K(t+) = P(OQ(t +) — Poet + Do! (+ DP() 


Leça Pope 


Kr+1) = POOCADÍ HOT + DPOt+D|- POA Do! (+ DPADOL+ 
Irqi +DP(Dq(t+ 1) 


Ka = PODE + O C++) q" (t+ POD) 
Lo! (t+ DPCOq(t+ 1) 


Eta) = ideia 
rot DPOp( +) 
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de 


identificação 


saída z 
sa —— 


UNO LA 


ka 


controle e 


o pe 


o) 100 200 300 


parâmetros estimados 


(o) 100 200 300 


(amostras) 


Exemplo 5.6 - Considere o processo discreto de primeira ordem 
nº 
Iá YO tay(t=1)=bou(t—1) +e(t) 


pn -0.8 e by = 1. À variância do ruído e(t) é 0.01. Admite-se que o 
E entrada é um sinal do tipo PRBS (Pseudo Random Binary Sequen- 
io" amplitude 7, Os dados de saída e entrada são ilustrados na 
tivo ga As equações do estimador dos mínimos quadrados recur- 
são filed (5.25), são utilizadas para estimar a; e bo. As estimativas 
tivas e na figura (5.4) quando 0(0)=0 e P(0)=10lxs. As estima- 
(EBORE algumas observações, convergem para os valores verdadeiros 
G etal. 2004; ASTRÔM; WITTENMARK, 2011). 
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Observações: 

= Avariável t+ 1) no 
e(t+1) - y(t+ 1) 
e a saída do processo é à predição da saída util 
y de parâmetros (1) a partir dos 
te de tempo t. Por essa razão, é 
de predição da saída ou série de 


algoritmo dos MOR 
q! (t+ DÃO 


é o erro entr 
zando as estimativas do veto 
instantes anteriores ao instan 
usualmente denominado erro 
inovações; a 

= O erro de predição torna-se zero quando 0(t) > 0. Assim, a 
medida da qualidade do estimador reflete-se na matriz defini. 
da positiva, P(t), denominada matriz de covariância do estima- 
dor À . O tamanho (magnitude) dos elementos da diagonal de 
P(t) está relacionado com a variabilidade dos correspondentes 
elementos em À. Por exemplo, se o elemento (1,1) é pequeno, 
significa que a estimativa de à, é adequada (variância baixa). 
Por outro lado, se o elemento (na + 1, na + 1) é grande, signi- 
fica que a estimativa é inadequada (variância alta); 


* Para a inicialização do algoritmo dos MQR, deve-se atribuir os 
valores de 6(0) e P(0). Se uma condição inicial para os valores 
dos parâmetros está disponível a partir do conhecimento à 


priori do processo, então estes valores devem ser utilizados 
para 0(0) e P(0) 
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» | processo 


I E 


entrada | | saida 


modelo 
estimado 


A cada período de amostragem, utilizam-se as medidas de en- 
trada e saída do processo e existe limitação no tempo de pro- 
cessamento; 


No contexto de controle adaptativo durante a fase inicial 
de processamento, quando os parâmetros do processo estão 
sendo sintonizados, o sinal de controle é inadequado e uma 
mudança é produzida na variável do processo, comprometen- 
do o comportamento dinâmico do sistema de malha fechada; 


Se as estimativas iniciais são pobres, então P(t) é construída 
como uma matriz diagonal com elementos positivos e de 
magnitude elevada. Quando as estimativas melhoram, os 
elementos de P(t) decrescem em magnitude, de modo que 
o ganho, K(t), torna-se aproximadamente nulo, resultando 
04 +1)= O): 

Uma vez parametrizado o processo, deve-se qualificar o mo- 
delo estimado utilizando técnicas de validação de modelos. 
Entre as diversas técnicas de validação, pode-se investigar a 
magnitude de certos índices de desempenho (BROSILOW; 
JOSEPH, 2002). Os índices de desempenho para a avaliação da 
qualidade dos modelos matemáticos podem ser calculados 
pelas seguintes equações: 
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E - 2 
seQ => yo) 49] 


Je 


Se é avaliado o somatório do erro quadrático médio, então 
SEQ/N. O modelo linear estimado que melhor se ajuste as me- 
didas gera o menor SEQ; 


» a att 2 
Coeficiente de Correlação Múltipla - R 


N 
> lydo = 309] 
R=I-Ho 
Sly) 
[e 
onde y(k) é a saída real, P(k) é a saída estimada e Y é a média 
das N amostras da experimentação. Quando o valor de Rº é 
igual a 1 (um), indica uma exata adequação do modelo para os 
dados medidos do processo. O valor de Rº entre 0.8 e 1 (um) 
pode ser considerado suficiente para muitas aplicações prát- 
cas em identificação; 
* Segundo Brosilow e Joseph (2002), tem-se uma adequada 
excitação de entrada na planta quando a razão sinal-ruído 
(S/R) é aproximadamente 5. Para obter uma estimativa d2 
variância do ruído, coletam-se N dados da planta em um i- 
tervalo de interesse. A razão S/R da saída da planta é 
pela seguinte equação: 


7) Veda) 


A 


onde &'/ é à variância da saída durante o teste de identificação é 
9; repre-sentaa variância da tante, 
sendo calculadas por saída para uma entrada cons 
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»y(t)-3] 
g)=g'=tl 


N 


. À qualidade do modelo estimado de um sistema depende em 
parte da natureza do sinal de entrada aplicado durante 4 fa- 
se de coleta das amostras. Isso acontece porque o sinal de 
entrada excita, em variados graus de intensidade que depen- 
dem de suas características, as denominadas frequências 
naturais ou modos do sistema, ou seja, o sinal de entrada força 
o sistema a revelar, na saída, as suas características dinâmi- 
cas. Assim, o sinal PRBS, ao ser aplicado como entrada no 
processo real, tem dupla finalidade: a) excitar os modos do 
sistema que correspondem ao conteúdo espectral; b) prevenir 
a ocorrência de um mau condicionamento numérico e/ou 
colunas iguais na matriz de covariância. O sinal PRBS só pode 
assumir dois valores (ou níveis), Zó. A comutação de um nível 
para outro ocorre de forma aleatória, cuja estatística pode ser 
controlada (seleção do conteúdo espectral do sinal PRBS) e 
baseia-se em determinadas regras de implementação (LJUNG, 
1999). Por exemplo, um sinal PRBS pode ser gerado a partir 
de um registrador de deslocamento, portas lógicas Pci 
do tipo E e OU, cujo comprimento máximo da agree E 
N, = 2% — 1, onde n, é o número de bits do registrador : 
locamento (BARREIROS, 1989). Um Cape cores Deise 
de implementação de uma sequência P RES rar uma se- 
lógica: utilizar a função rand do Matlab para e 
quência de números aleatórios uniformement Definir uma 

intervalo [0, 1), armazenando-a no vetor mesh 

variável alfa com um valor qualquer entre 

cada elemento do vetor sna maior ou mpi outro vetor de- 
valor +ô ao elemento correspondente — vetor sna menor que 
nominado PRBS; para cada elemento a te do vetor 
alfa, associa-se o valor -Sao elemento sna são uniformemente 

PRBS. Como os elementos do vetor 
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distribuídos no intervalo [0, 1), o valor médio dos elemento, ; 
0.5. Portanto, uma seleção de alfa maior que 0.5 Proporciona 
um nível tendencioso para -ó, enquanto que a seleção de alfa 
menor que 0.5 proporciona um nível tendencioso para , ê 
Mediante a seleção de alfa, pode-se induzir uma maior a 
menor taxa de comutação entre +ôe -ó, controlando assim o 
conteúdo espectral do sinal PRBS. Aumenta-se a frequência 
da sequência PRBS (taxa de comutação entre +0e 9) se 
valor de alfa está próximo de 0.5. A amplitude do sinal PRBS é 
alterada através de um fator multiplicativo beta (SAKAI, 1999) 
A tabela (5.4) implementa em Matlab a sequência PRBS para 


alfa= 0.5, beta= 1etf= 5. | 
Tabela 5.4 - Programação em Matlab de uma PRBS ] 
$ Sequência binária pseudoaleatória t 
clear all 
delta=input ('Intervalo de integração DELTA t = '); 8 
tf=input('Tempo final tf ="); 7 
t=[0:delta:tf];sna=rand(size(t));1=0; 7 
while 1==0 À 
alfa-=input('Valor de alfa = '); h 
if alfa<=1 q 
if alfa>=0 
1=1; 
end 
end 
if 1== 
disp('Alfa deve estar entre Q e 1'); 
end 
end 


for i=l;length(t), 
1f sna(ij<=alfa, sbpa(i)=-1; 
else sbpa (i)=1; 
end 

end 

plot (t, sbpa); 


dócitt bra a A 


Outra sugestão de projeto do sinal de entrada do tipo PR$S 
conforme discutido em Rivera e Flores (2000) e Brosi 
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Joseph (2002) para sistemas mMonovariáv: z 
Inicialmente o usuário deve ter uma estimativa d 
tes de tempo rápida (menor), 1, e lenta (maior) xe Tapaç 
terizam o sistema sob análise. Assim, define-se & 
banda de frequências para excitação do sinal de en 
também está associada à parte plana do espectro s 
o sinal de entrada tem potência suficiente), isto é, 


eis, é apresentada. 


que carac- 
largura de 
trada (que 
obre o qual 


onde /% é o fator representando o tempo de acomodação do 
processo (para Tosz, /k= 3; para Tox, = 5) e q é o fator 
representando a velocidade da resposta de malha fechada 
como um múltiplo do tempo de resposta de malha aberta. A 
geração do sinal PRBS é caracterizada por dois parâmetros: o 
número de bits do registrador, n,, e o tempo de chaveamento, 
Tu. À sequência repete-se após N.T., unidades de tempo, onde 
N.=2"— 1. A caracterização do tempo de chaveamento e o 


comprimento da sequência obedecem às seguintes relações: 


24 
E: < 2.87, E N =92% E cuia 
A Io 


onde n, e N, são valores inteiros e Ts deve ser um pena 
inteiro do período de amostragem T,. Admitindo que % = e 
min, /, = 5, o, = 2, obtêm-se Tw = 14, Mr =& NsTow a 
Para uma amplitude de 42.5, o seguinte comando bai: 


pode ser executado (BROSILOW; JOSEPH, 2002): u = 


9, T, amplitude) = prbs (0.02, 0.2, 1.4, 2.5); he 
Uma excitação consistente em identificação de nr de excitação 
que o sinal de entrada satisfaça 0 regue E tação estatis 
Persistente, Adicionalmente, à PES dep oferece 
tica entre os sinais de entrada, saída é 
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algumas vantagens práticas o ea Significa 
procurar o melhor intervalo de me idas prai entáficação da 
planta, observando se uma das seguintes an isa residuais é 
satisfeita: a) calcular a função correlação cruzada entre os 

sinais de entrada e erro de predição pela seguinte equação; 


rud) = Setou(t+ k) 


O modelo estimado é adequado (os sinais de entrada e erro de 
predição não são correlacionados) quando os valores de r.(k) 
parak=0, 41, £2, £3, ... são aproximadamente nulos. Em 
aplicações práticas, utiliza-se um intervalo de confiança de + 
10%; b) calcular a função correlação para o erro de predição da 
saída pela seguinte equação: 


r(k)= Setvett +k) 


As medidas são adequadas para modelagem quando o conjun- 
to de elementos de r«k) para k = O está no intervalo de confi- 
ança de + 10% do valor normalizado r(0) = 1 (ISERMANN, 
1980; BROSILOW; JOSEPH, 2002); 


Se os parâmetros são variantes no tempo, deve-se modificar o 
estimador para evitar que o ganho torne-se pequeno, manter- 


do-se assim a capacidade de adaptação do estimador dos MOR 
(estimador em alerta). 


5.5 ESTIMAÇÃO DE PROCESSOS VARIANTES 
NO TEMPO 


O rastreamento de parâmetros variante é um impor” 
tante problema na identificação de sistemas dinâmicas, Independer” 
temente da aplicação dada ao estimador, é necessário que ele seja sie 
no rastreamento das mudanças no sistema identificado (presença de não 


ENADE Se 
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put — 
idades). Adicionalmente, deve-se obter um modelo 
poe resente adequadamente o comportamento 


Matemático 
do sistema coi 


que rc à aplicação em controladores adaptativos (WELLSTE ntrolado, 
a CLUETT, 1991). LISTED; ZARROP, 
199%; 


« Como aumentar a adaptabilidade do algoritmo dos MOR? 


5.5.1 PROBLEMA DO RASTREAMENTO DE 
PARÂMETROS VARIANTES 


Quando o número de iterações aumenta, os parâmetros estimados 
convergir. Esta convergência é normalmente refletida pela di- 
minuição dos elementos da matriz P(t) e é desejável para o caso de sis- 
temas invariantes no tempo. Entretanto, quando o sistema é variante 
notempo, é necessário evitar que os elementos da matriz P(t) se tornem 
“muito” pequenos, para que possa ocorrer a correção dos estimado- 
res 1). Na prática, procura-se um compromisso entre a capacidade de 
adaptação (P(t) grande) e a convergência no algoritmo de estimação 
(P(t) pequeno). 
Para manter a capacidade de adaptação do estimador, limitando a 
valores mínimos os elementos da matriz de covariância (controlar a 
Zagnitude dos elementos de P(t)), dois procedimentos são avaliados: 


= Atualização de P(t) e 
= Fator de esquecimento. 


3.5.2 ATUALIZAÇÃO DE P(t) 


: P(o) 
Para a atualização da matriz de covariância, os elementos de "o 
Fº aumentados pela adição de uma matriz diagonal semidefânida por 


“ad, Então, a nova matriz de covariância é calculada por 


P(D=P(D)+Q() 


tos basea- 
dog Aliteratura de identificação apresenta dois pre 
(5.26), isto é, random walk e covaria! 


(5.26) 


[a 
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TE INÂMicos 
5.5.2.1 BUSCA ALEATÓRIA [RANDOM WALK) 


ca random walk, adiciona-se, a cada iteração 


Para a aplicação da técni icásd 
N o 5 os' ZãO 1 
a matriz Q(t) cujos elementos descrevem a sup! € Variação dos 


parâmetros do sistema. Por questão de simplic idade e estabilidade do 
algoritmo dos MQR, usualmente a matriz Q(t) apresenta estrutura dia 
gonal, Q(t) = qUnasnba1)x(narmbel) onde q e *. a 

Se somente certos parâmetros mudam, então O(t) pode ser seja. 
cionada com zeros na diagonal em todas as posições exceto naquelas que 
correspondem aos parâmetros variantes. Por exemplo, se é conhecido 
que somente o primeiro e O terceiro elementos no vetor de parâmetros 
modificam-se no sistema com quatro parâmetros, então Q(t) é dada por 


Q(t)=diag(q;, 0, 92,0) 


sc ra sa 


Os valores de q; e q: devem refletir a magnitude da mudança dos 
parâmetros. Por exemplo, se o parâmetro 8, muda “pouco”, então q; = 0.05 
deve ser tentado; se o parâmetro 6: muda em 100%, então qs = 0.5 deve 
ser tentado. Além disso, existe considerável liberdade na seleção do 
tamanho dos elementos da diagonal de Q(t). Deve-se lembrar que são 
positivos porque P(t) é uma matriz definida positiva. Também é 'm 
portante observar que Q(t) aumenta o comprimento do passo de ajus* 
e um valor “grande” causa uma excessiva oscilação nos correspondent** 
parâmetros estimados. 


A equação (5.27) foi proposta e testada heuristicamente por Aruês 
pg e pode ser utilizada na seleção da matriz Q(t) = aUnaenbe tnaemdelh 
isto é, 

qe AO] e [200] 61 
[na+nb+1]) p 


onde t,[P(t)] é o traço da matriz ont 
mero de parâmetros estimados. P(t) ep = (na + nb + 1) representa ; 


A influência da técnica rando; amét” 
m walk nas estimativas par 
cas pode ser observada no exemplo 5.7 Pinda a EA OP, 199) 
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gaemplo 5.7- Considere um sistema caracterizado pelo seguint do 

e modelo; 
E y(O=-açy(t= D+ bçu(t=1)+b u(t=2) re(t) 
dos [0.7 se 0 g 

+ se 0<t<200 
=1,b=0.5e a 

ad onde bu = 1, 1 * 0.35 se t>200 


A figura (5. 6) ilustra os parâmetros estimados utilizando-se 6 
ele. NOR: «padrão com 9(0)=0, P(O) = 10013 e sem a aplicação do proce 
dimento random walk. Na figura (5.7), o método random walk é aplica- 


ad, 
[em do no instante de tempo t = 200 quando a; varia, onde O(t) = 10kss e 
por Qt) = 0.1Ias, respectivamente. 

Figura 5.6 - Estimador dos MOR sem a técnica de busca aleatória 

1.5 — E as - | 
tdos a 
0.05 JE eme 
deve 
pus 05 
2 são 
Fim 6 
juste 9 
entes 

-05 
quiz 
mbelh EA 

É aoo 
5.20 aço 2007-8004 0400: 08. 609 
a amostras 

se 
pó 6) e (5.7), observar 
á à De acordo com os resultados das figuras di ! uma 
get à utilização do estimador dos MQR a a eg 
390» “pveEência lenta nos parâmetros estima ações nas estimativas nd 
” para Q(t) conduz a excessivas variaç ra Q(t), tem-se 


pa 
Fr que, utilizando-se uma matriz “pequena a aplicação da t 
transitório nos parâmetros estim: 
a 
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DS Da 


random walk deve-se conhecer a priori, para à aplicação particular, que a, 


é o parâmetro variante do processo. 
atória: [a] Q(t) = 101355, [b] QUO) = 0.11s5 


MQR com busca ale 


F 
gur 
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5.5.2.2 REINICIALIZAÇÃO DE P/t/ICOVARIANCE RESETTINGI 


Se a matriz Q(O) é adicionada com muita frequência, os elemen 
tos da matriz P(t) tornam-se grandes, causando uma flutuação nas 
os da matriz 
etimativas. Por outro lado, se a frequência é pequena, os elementos da 
estimativas 


matriz P(t) tornam-se pequenos à ponto de o algoritmo não rastrear as 
mudanças no sistema (o algoritmo “adormeceu”). Portanto, existe um 
instante (intervalo) de tempo ótimo para a alteração de P(t). A técnica 
de reinicialização da matriz de covariância evita este problema pela 
adição da matriz Q(t) somente se a magnitude dos elementos de P(t) é 
menor que um determinado valor especificado pelo usuário. A magni- 
tude da matriz de covariância é avaliada ou pelo traço de P(t) ou pelo 
erro de predição. 

Dado um traço mínimo to, nas iterações onde t,[P(t)] « tro, faz-se 
uma reinicialização de P(t) somando a matriz Q(t) cujos elementos dão 
um ganho abrupto ao estimador. 

Outro procedimento para a regulagem dos elementos de P(t) uti- 
liza o erro de estimação /v(t)- 5(t)], que também pode ser observado e 
implementado para a reinicialização de P(t) somente se o valor absolu- 
to do erro de estimação excede um valor especificado pelo usuário. 


5.5.3 FATOR DE ESQUECIMENTO 


Para processos variantes no tempo, é necessário fornecer ao 
algoritmo dos mínimos quadrados recursivo uma capacidade de adap- 
“ação mínima, impedindo que o ganho do estimador tenda a zero. Esta 
capacidade pode ser obtida dando-se uma maior importância às novas 
medidas pela regulagem do fator de esquecimento. 


*% Introduzir uma constante no algoritmo (denominada fator de 
esquecimento) que pondera mais as últimas medidas. 


O algoritmo dos MQR-padrão foi desenvolvido considerando-se 


e. Plt+ 1) = [hp (te D(t+ 1) 


He Dt 1) = "(OM + (t+ DO! (t+ 1) (5.28) 
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SINE 
Se os parâmetros são variantes, uma modificação no esti 
dos MQR faz-se necessária para evitar que o ganho torne-se * 


pequeno. Assim, em vez da equação (5.28), pode-se utilizar 


dra Dto = AO! COPO PCt+ DO! (t+) 


Mador 
Muitos 


(5.29) 
onde À é o fator de esquecimento, assumindo valores entre () (zero). 
1 (um) e penalizando as informações antigas na presença de novas mp. 
didas. 

O algoritmo dos MQR-padrão minimiza o seguinte critério: 


N 
IN9)=5 (0-9 (00]º 


t=1 


em cada instante de tempo e pondera cada erro de estimação, e(i) = yíi)- 
((1)8, igualmente. Uma mudança paramétrica no sistema identificado 
significa que erros recentes devem ser mais fortemente enfatizados do 
que erros anteriores de modo a permitir ao estimador dos MG 
adaptar-se ao novo sistema. Assim, para modelar as variações dos 
parâmetros variantes, o critério dos MQR modifica-se para 


J(N,9)= Deo 


t=1 


e o efeito do desconto dos erros passados pode ser visto como 


J(N,0)=2JK(N-1,0)+e?(N) 


Deste modo, as medidas velhas são i “esqueo” 

E exponencialmente je 
das” e maior ênfase é atribuída às novas medidas. O algoritmo das 
mação dos mínimos quadrados com fator de esquecimento apresen? ; 
seguinte forma recursiva (LJUNG, 1999): 


Der = AO ++) qr) = q" (t+ DÃO) 


KtsD=—  POgt+) 


E (5.30) 
A+ qt DPOp(t+1) 
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| 1) P(tp(t + Dp"! (t+ PC) | 
, == p(y) —— ZA es 
O CEDRO) 


Observação: 
* Oalgoritmo dos MOR com fator de esquecimento evita que os 
Ê elementos de P(t) tendam a zero, mantendo o estimador em 
: “alerta” para rastrear dinâmicas variantes, isto é, 
se À=1, tem-se a mesma ponderação (MQR-padrão); 
— na prática, tem-se 0.9 <2 <1 (ponderação maior para as me- 
didas atuais). 


Para aumentar a sensibilidade do estimador dos MQR na presença 
de variações nos parâmetros do sistema, implementa-se um fator de 
esquecimento variável. Esta técnica pode utilizar o erro de predição para 
verificar, em cada período de amostragem, o estado do estimador. Se o 
erro é menor que um determinado valor fornecido pelo operador, então 
os parâmetros estimados estão próximos dos valores verdadeiros. Neste 
caso, o fator de esquecimento deve assumir valor próximo de um. Se o 
erro é grande, o fator de esquecimento deve assumir valores pequenos 
(0.9<2 <1) de forma que a sensibilidade do estimador aumente para 
permitir o rastreamento dos parâmetros e, consequentemente, reduzir 
o erro de predição. Entre os vários métodos para a regulação do fator de 
esquecimento, tem-se 


al Fator de esquecimento constante 
Neste método, A(t) é calculado por 


| A(t)= À, onde 0.9SA, <1; 


b) Fator de esquecimento variável 
Fator de esquecimento com traço limitado de P(t) 


Neste método, calcula-se A(t) para manter o traço da matriz de co- 
“ariância P(t) superior a uma constante ty definida pelo usuário, isto é, 
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| e feco] , se efero to 
Mt)= I 


to 
1, de outra forma 


Seleciona-se O traço mínimo (to) de modo que, em regime, o = 

dor apresente pequena variância em torno do valor estimado, e 

da deve permitir uma adaptação a variações lentas nos par a 
tros a serem identificados; 


Fator de esquecimento exponencial 

Neste método, procura-se manter At) em 1 ou na vizinhança 
(Alt) = 0.999). Quando uma mudança paramétrica é detectada, i(t) 
pode ser ajustado para 0.95, o que implica uma rápida variação nos 
parâmetros estimados, e o comportamento de Alt) tem um ajuste expo- 
nencial crescendo para 1 ou 0.999. Assim, é possível ajustar 2(t) de 
acordo com 


MO= MAMA A 


onde 25 < 1 e 2(0) < 1. Por exemplo, para 2, = 0.95 e (0) = 0.95, resulta 
que 2(5) = 0.9632; 1(10) = 0.9715; (20) = 0.9829 


Mb= Aoh(t-1) +. [1 —Ao] 
N 
este caso, para Ày = 0, 95eÃ= 0.999, obtém-se 


M=0.95A(t- 
&, em regime, tem-se Mt-1)+0.999(0.05) 
limA(t) = = 0.999 


Digitalizado com CamScanner 


07 25 


ados por equações a diferenças 161 


s Identificação de sis 


A influência da técnica do fator de esquecimento nas estimati- 
aramétricas pode ser observada no exemplo 5.8 (WELLSTEAD: 
as 


ZARROP, 1991). 


Exemplo 5.8 — Considere um sistema paramétrico representado por 
y(D=-a,y(t-1)-asy(t-2)+bou(t-1) +e(t) 


As figuras (5.8) e (5.9) mostram a influência sobre o estimador dos 
MGQR para diferentes valores do fator de esquecimento. Neste exemplo, a 
o coeficiente by muda em t = 200 e os coeficiente a; e a; são constantes, 
isto é, 


a=-1; a,=0.25; b=1(t<200) ; b,=2(t>200) 


Observe que o efeito de redução de A(t) aumenta a sensibilidade 
das estimativas de mudança do parâmetro bo. Uma penalidade associa- 
da, entretanto, é que a variabilidade das estimativas aumenta. 


Figura 5.8 - Estimador dos MQR com 2/t) = 0.99 
25 


o 100 200 300 400 500 600 700 800 
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Figura 5.9 = Estimador dos MaR com: la) Aft! = 0.98. (b) Alt) = 097 


25 


o 100 200 300 400 500 600 700 80 
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Observações: 

* Autilização do fator de esquecimento variante no tempo per- 
mite que o algoritmo dos MQR rastreie melhor as variações 
nos parâmetros do que o algoritmo com fator de esqueci- 
mento constante. Entretanto, com um fator de esquecimen- 
to constante ou variável, aumentam-se de um mesmo valor 
todos os elementos da matriz P(t). Assim, os parâmetros que 
não são variantes no tempo oscilam; 


= Através da escolha apropriada dos elementos da matriz Q(t), 
na aplicação da técnica da atualização de P(t), os elementos da 
matriz de covariância e o ganho do estimador podem não 
tender mais a zero. Este método é mais seletivo que a utiliza- 
ção do fator de esquecimento, pois permite aumentar separa- 
damente os elementos individuais da matriz P(t); 


* Para um dado fator de esquecimento, o estimador de parâme- 
tros conduz a estimativas baseadas nas últimas N2 amostras 
de dados (memória do estimador), onde 


N, = as (constante de tempo associada com À) 
Assim, o valor de A(t) = 2o = 0.98 conduz à estimativa de pará- 
metros baseada essencialmente nas últimas N; = 50 amostras 
de dados, enquanto que o valor de o = 0.995 conduz à 
estimativa de parâmetros baseada essencialmente nas últimas 
N; = 200 amostras; 


* O procedimento do fator de esquecimento funciona de forma 
adequada somente se o processo tem excitação. O uso ina- 
dequado do fator de esquecimento às vezes produz uma 
condição inadequada, denominada blow up, na matriz de 
covariância. Isso acontece quando um fator de esquecimento 
constante A(t) « 1 é utilizado e o sistema está operando em 
regime por um longo período de tempo (pode também acon- 
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5, 0). Neste caso, a matriz 
tecer se P(t As ! seguint E AUment; 
exponenc jalmente, de acordo com à seguinte equação 
P(t) P(t-1) 
(t) = 
Mt) 


e evita-se este problema de operação do estimador dos Mgp 
com um fator de esquecimento variável. A figura (5.10) ostra 
a presença do fenômeno blow up a partir do instante 409 
quando o traço da matriz de covariância cresce significa 


vamente; 


ESA grrsa va: S 


Figura 5.10 - Presença de blow up pela falta de excitação 
convergência dos parâmetros estimados 


5 


ESSE 


Fs 


é Rd 
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A utilização de algoritmos de estimação paramétricos em com- 
putadores, para ser bem-sucedida, exige que o algoritmo seja 
robusto o suficiente para contornar os problemas que fre- 


quentemente ocorrem nas implementações em tempo real 
(robustez numérica do estimador). Os problemas relacionados 
às limitações dos equipamentos (conversores, sensores, erros 
de arredondamento dos computadores) introduzem impreci- 
são numérica devido ao mau condicionamento da matriz de 
covariância que perde sua propriedade de matriz definida 
positiva, instabilizando o algoritmo de identificação. Este pro- 
blema é solucionado se a matriz de covariância é decomposta 
de acordo com 


P(t)=UCOD(YU” (t) 


onde U(t) é uma matriz triangular superior com todos os seus 
elementos da diagonal unitária e D(t) é uma matriz diagonal. 
Na literatura, este método é conhecido como fatoração UD 
(Upper Diagonal) de Bierman (1977) e garante que a matriz 
P(t) seja definida positiva e que tenha boa estabilidade numé- 
rica com erros numéricos diminuídos. Ástrôm e Wittenmark 
(1995) apresentam o código em Pascal do estimador de pará- 
metros dos mínimos quadrados utilizando a fatoração UD; 


Para solucionar os problemas de estimação de parâmetros, é 
necessário inserir no sistema um módulo de supervisão basea- 
do em procedimentos heurísticos. Desta forma, funções de 
segurança são implementadas no nível superior para melho- 
rar o desempenho dos algoritmos de estimação no nível 
inferior, a partir das medidas do processo sob identificação 
(figura 5.11). Em processos complexos, a presença de parâme- 
tros e perturbações não modeladas pode comprometer o 
desempenho dos algoritmos de identificação. Para prevenir 
imprecisões no algoritmo de estimação com respeito a tum 
off, windup, parâmetros variantes no tempo e dinâmicas não 
modeladas, o módulo de supervisão deve ser implementado 
para: a) monitorar o erro de predição e estimativas de pará- 
metros para mudanças abruptas no processo que resultam na 


[a 
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ruptura do modelo estimado; b) ajustar as medidas e os EE E 
metros de projeto do algoritmo de estimação a mudanças à 
condições de operação da planta (ISERMANN, 1980); Ty 


Figura 5.11 = Estimação com supervisão 


módulo de 
supervisão: 
— monitoração 
— diagnóstico 


saída 


, A determinação da ordem do sistema é uma importante tarefa 
na estimação de parâmetros (HANG et al., 1993). Na apresen 
ne dos exemplos anteriores, foi assumido o conhecimento 

ordem dos processos, isto é, o valor de na = n era conhecido 

qe prática, entretanto, a ordem do processo pode não 

pe em conhecida. Se um modelo de ordem incortetê 

pio pd ri ocorrer erros ou redundâncias de termo* 
peida mpo de processamento do algoritmo). 

Eng modelo empregado na estimação de parâmetros é 
+ à soma dos quadrados do erro de predição 


Cho Hbo-e! pô) 


diminui à medida 
né utilizado ia Que se obtém o melhor ajuste. O subsar 
Para denotar o valor do e 


Digitalizado com CamScanner 


os Identificação de 


sistemas representados por equações a diferenças 167 


ordem n previamente selecionada. Entretanto, se a diminui- 
ção de S, é lenta entre quaisquer dois modelos, com o segundo 
modelo de ordem maior que o primeiro, então a utilização de 
um modelo de ordem maior não reduzirá significativamente 
S.. Assim, para determinar a ordem do sistema de um 
processo, certo conjunto de medidas, isto é, um conjunto das 
variáveis da planta (u(t)) e fy(t)] pode ser utilizado para a 
estimação de parâmetros usando um modelo de ordem 
particular. O mesmo procedimento pode então ser conduzido 
com os dados utilizando modelos de ordens diferentes. 
Observando o gráfico de S, com a ordem do sistema n, uma 
ordem do modelo pode ser selecionada, por exemplo n', onde 
a inclinação 48,/4n é íngreme para n < n' e suave paran > nº. 
Este procedimento está ilustrado na figura (5.12); 


Figura 5.12 - Teste de detecção da ordem do sistema 


O atraso de transporte d do processo pode também ser iden- 
tificado de forma similar ao caso da determinação da ordem 
do modelo. Para uma dada ordem do modelo e uma sequência 
de valores de d (por exemplo, d = 1, 2, ...), a melhor estimativa 
de d é a que conduz ao menor valor da soma dos quadrados do 


Bórtme 
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erro de predição. Portanto, é possível combinar a determina. 
ção do atraso de transporte e a ordem do modelo do Processo 
em um problema de identificação pela simples interação dk 
procedimento de detecção do atraso para uma sequência de n 
valores (HANG et al., 1993; MORAES, 2004); 


* Critérios de Informação (CI) podem ser avaliados para manter 
a ordem do modelo da planta estimada tão simples quanto 
possível (princípio da parcimônia). Critérios como de informa- 
ção Bayesiana, de Akaike ou minimum description length com- 
binam a variância residual do erro de estimação e a ordem do 
modelo (HABER; UNBEHAUEN, 1990). A figura (5.13) ilustra 
o comportamento do critério de informação para detectar a 
ordem do modelo estimado com base no menor resíduo de es- 
timação associado ao menor conjunto de parâmetros estimados. 


ssa 


e 


Figura 5.13 - Evolução paramétrica do critério de informação 


critério de 
informação 


M, M, M, M modelo 


seres um modelo de baixa ordem é selecionado, 
exemplo ordem 1 (um), então aumenta-se a ordem do 
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apiT 


mado e o CI é avaliado para cada incremento na ordem do 


estir 
modelo. Assim, a estrutura adequada do modelo estimado é a 


f. 


que proporciona a menor taxa de variação do critério de infor- 
mação. Diversos critérios de informação podem ser obtidos da 
literatura de identificação como uma complementação da 
função custo básica dos mínimos quadrados, isto é, 


É 


FÊ 


N 
Ju =D yd9- Pd 
N& 


Ma, 
Mem : ; 
com alguns termos extras que penalizam a complexidade do 
h da modelo matemático da planta e podem ser selecionados de 
acordo com um dos seguintes testes (HABER; UNBEHAUEN, 
[e 1990; LIUNG, 1996): 
im - Critério de Informação de Akaike (Akaike Information Criterion) 


Ea 


AIC=Nin[3,]+2p ou ac=(14 833, 


onde N é o número de medidas da experimentação e p é o 
número de parâmetros do modelo estimado; 


— Erro de Predição Final (Final Prediction Error) 
N+P! 4 ppp=(N+P) 


N-p] (N-p) 


EPE=Nin(3,J+ Na Fe 


— Critério de Informação Bayesiana (Bayesian Information Criterion) 
BIC=NIn[J,]+pln[N]. 


A 36 ALGORITMO DE ESTIMAÇÃO DA 
APROXIMAÇÃO ESTOCÁSTICA 


É | x O estimador dos mínimos quadrados recursivo apresenta um 
A | Úmero significativo de cálculos a cada iteração, o que pode inviabiizá-lo 
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em algumas aplicações práticas envolvendo controle adaptativo a 
patibilizar o tempo de processamento em tempo real do algoritmo E 

a dinâmica da planta a ser controlada). Isso se deve ao fato de o Processo 
apresentar uma constante de tempo pequena (processo rápido) ou - 
processo é complexo, isto é, se as ordens de na (associado aos Polos de 
malha aberta) e nb (associado aos zeros de malha aberta) são Brandes 
Assim, na equação de atualização das estimativas à 


Ot 1) = 0) + P(A Da(t+ 1e(t+ 1) (5a) 


deve-se substituir P(t + 1) por uma expressão que seja mais simples de 
calcular. Portanto, pode-se utilizar P(t + 1) na equação (5.31) pelo ajusta 


de um ganho xt + 1), ou seja, 
t+) = KO + (t+ Do(t+De(t+ 1) (530) 
cuja seleção obedece a um dos seguintes procedimentos: 
a) Sequência ganho escalar constante 


Y()=7, =escalar constante 


b] Seguência ganho com aprendizado sequencial 


(t)= RETO gra 
Bro" (ot) 


onde 20 e 0< «2 são escalares fixos ajustados pelo usuário pat 
calibrar a e: i 1º tamanhão do passo de adaptação do estin” 
qe : tificação de Sistemas Dinâmicos, este esti 
denomina-se algoritmo da Projeção normalizado. 


“| Sequência ganho com o vetor total de medidas 


t be 
"o-[Bao '| =[(0]" 68 
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Com a seleção de xt) pela equação (5.33), o algoritmo da equação 


(5.32) toma-se 


d0=0t=D + OO q" DO] ; 
r(bD)=r(t=—1) + loco 


Observações: 


» As vantagens do método de estimação da aproximação esto- 
cástica são a simplicidade e a velocidade de cálculo. As desvan- 
tagens do método de estimação da aproximação estocástica 
são a baixa convergência dos parâmetros estimados e a “pobre” 
adaptação a mudanças nos parâmetros variantes do sistema 
(WELLSTEAD; ZARROP, 1991); 


* No caso da sequência ganho escalar constante o algoritmo de- 
nomina-se LMS (Least Mean Squares) (LJUNG, 1999); 


* Outros algoritmos para identificação de sistemas dinâmicos 


utilizando a técnica da aproximação estocástica são (SINHA; 
KUSZTA, 1983): 


— Saridis e Stein (1968): 
h(o = ht =D + Du ID -u”(Dh(t=1)] 
onde 


h=[h,h,.ch]" 
u(o) =[u(t=1)u(t=2)... u(t=0F' 
YD=1/t(t=1,2,.0) 


e sendo h(t) a estimativa dos elementos da resposta impul- 
siva. Se a função de transferência discreta do modelo do sis- 
tema contém (na+nb+1) parâmetros, então /=na+nb+1 
amostras da resposta ao impulso são estimadas; 
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Kwatny (1972) É 

" ' 15 (y Le cen D)= y(OKt=1) 

n CE ms 
edi? jóce-2) 


onde a notação é similar ao estimador dos mínimos Quadra, 
dos e é uma constante positiva; 

» Se e(t) é um ruído branco, o estimador dos MQ é assintori 
camente não polarizado. Entretanto, no caso onde e(t) não é 
branco, tem-se a seguinte representação para a planta (mode. 
lo geral): 


Az y(t)=Blz Dult-d) +e(t) ; eD=C(z"w(t) 


onde w(t) é uma sequência ruído branco (Teorema da Represer- 
tação). Portanto, neste caso, deve-se verificar como alterar o 
estimador (A(z "),B(z"),C(z")) para que as estimativas sejam 
ainda não polarizadas (ASTRÔM, 1970). 


5.7 ALGORITMO DE ESTIMAÇÃO DA VARIÁVEL 
INSTRUMENTAL 


o Processo, e(t), é uma sequé 
os do estimados não polarizados e consistentes, 
trumental (VI) Petr é branca, utiliza-se o método da variável ins- 
denominada variá Por Young (1970), Assim, uma nova variá 

cioná-la com as dá "tal, é definida, procurando-se correlas 
ná-la do ruído do sistema, entrada e saída do processo e descorrelacio” 


vel instrume 
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A seguir, discute-se a aplicabilidade de implementação do método 

avi substituindo-se a equação (5.7) na equação (5.13), obtém-se o se- 
da VÊ. « 


e valor para a estimativa 0: 
0=[4"6] "o" [y0+E] 


Após algumas manipulações algébricas, tem-se 


d=0+[676]" [97] (5.34) 


A estimativa 6 é dita consistente se converge para O quando o 
número de amostras aumenta (N — q). Portanto, determina-se a 
consistência do estimador, equação (5.34), calculando a esperança ma- 
temática da variável À, =[/0], em termos de probabilidade, por 


=10)=0+=]676]"=[y"E] 


Para avaliar se Ô converge para 8, calculam-se os termos 5/6" 6] 
e=[6 E]. Se e(t) é uma perturbação do tipo ruído branco, o termo 
=/9 E] =0, a estimação é não polarizada e, consequentemente, ô con- 
verge para 0. Se a perturbação não é branca, o termo Z/9'E/+0, e o 
método dos mínimos quadrados fornece estimativas polarizadas, isto é, 
Í não converge para, inviabilizando sua implementação (GOODWIN; 
PAYNE, 1977). Este problema de polarização pode ser contornado com a 
técnica da variável instrumental. Assim, deve-se selecionar uma matriz 
Variável instrumental, Z, composta de elementos da variável instrumental 
(novo sinal), z(t), tal que z(t) seja correlacionada com as medidas do 
Processo (entrada e saída) e não correlacionada com o ruído, isto é, 


=lz'e)= 0; =lz"y|- W (matriz não singular) 


(alg À estimativa de , obtida pelo método da variável instrumental 
ritmo não recursivo), é calculada conforme a equação (5.35), ou seja, 


du=[2" 00]! Z' vm (5.35) 
Onde à matriz Z(t) é dada por 


guint 


[a 
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2"(1) 


Z(t) = 
2(N) ] 


| 


e dt) e Y(t) estão definidos nas seções anteriores. À análise da conyer. 
gência do algoritmo da variável instrumental está disponível em Ljun, 

(1999). Muitas vezes, à dificuldade em selecionar a variável instrumen. 
tal z(t) torna o método da VI não atrativo para determinadas aplicações 
(o algoritmo da VI não estima os parâmetros do modelo do ruído). 


5.7.1 ALGORITMO NÃO RECURSIVO DA VARIÁVEL 
INSTRUMENTAL 
Considere novamente o modelo da equação (5.7), isto é, 
Y=40,+E 


onde G representa os parâmetros verdadeiros do sistema que devem ser 
estimados. Multiplicando-se pela matriz da VI transposta, Z”, tem-se 


que 
Z'Y=2"40,+2Z7E 
ve =feofevl- fo) ee] 


Admitindo-se que a equação 


du) =[z"y]" [27] 

=0,+[2"9)"[z"2] 

calcula o algoritmo da VI não 

7 7) recursivo (mét ine), então à com 
dição 6, (N) > 0,, quando N -» 0x, é ass método off-line), en! 

hipóteses de que egurada se está subo 


O sinal 2(t) não está Sorrelacionado com e(t); 
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+ Osinale(t) tem média nula; 


+ Ainversa da matriz [z" 6] * existe, 


5.7.2 ALGORITMO RECURSIVO DA VARIÁ 
INSTRUMENTAL BHAVEL 


As equações para o estimador recursivo da variável ; 

em ser obtidas de maneira similar às e pç 

reursivo (LJUNG; SODERSTRÓM, 1983). A partir das Efe 

entrada e saída do processo, obtém-se o algoritmo da variável ng 
mental recursivo (método on-line) segundo as seguintes equações: 


Ke+1)= 500 +K(+D pt +) -q"(e+ não) 


P(t)z(t + 1) 
Kt+D)= às 
E raro) CDMA (629 
1) Pltz(t+ 1)” (t+ neo) 
P = = 
Sd y tádi A+ q! (t+ 1)P(Dz(t+1) 


Observações: 
* As mesmas considerações adotadas na inicialização de X0) e 


P(0) no estimador dos MQR podem também ser aplicadas no 
algoritmo da variável instrumental; 
de esquecimento como a 


. 
De modo análogo, tanto o fator 
atualização de P(t) podem ser implementados É a 
lhorar a capacidade adaptativa do estima 
Sistemas variantes no tempo. 


p AL 
5.7.3 SELEÇÃO DA VARIÁVEL INSTRUMENT tal 
o instrumental, 
“0, Diversos procedimentos para a seleção es ri proces- 
dos (sçm sido apresentados na literatura de e ementos de z(t) são 
DERSTRÓM et al., 1987). Geralmente, ae 
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los como Os valores atrasados ou filtrados dos sinais de entr 

a melhor compreensão, seja a representação do dE 

i ; ia- 

ão da VI, a partir do modelo nominal da pica 

a, 


selecionad 
da e saída. Para um 
grama de blocos da seleç 
apresentada na figura (5.14). 

1 Diagrama esquemático da variável instrumenta 


Figura 5 14 
e(t) 


modelo nominal 


u(t) 


A seguir, discutem-se algumas abordagens fundamentada: 
trabalhos de Young (1970), Sôderstró i Welstead é 
Zarrop (1991) e Ljung (1999). END E MN 


Considere a seguinte variável instrumental: 
t)=[-x(t =1)...-x(t-na) u(t -d)...u(t-d —nb)]" (537) 
onde o sinal x(t) é obtido do modelo auxiliar dado por 
At) =Blz Yu(t-d) o 


e os polinômios Z/: E 

e nb Ape e og ig po ") devem ser estáveis, de ordens nº 214 

variável instrumental de acordo tros reais do sistema. Seleciona-se à 
Srdo com um dos seguintes procedimentos: 


a) Variável instrumental via filtros adaptativos 


Quando os valo; 
então obtêm-se as E amétricos do sistema são desconhecidos; 
estimativas dos parâmetros do sistema 


Digitalizado com CamScanner 


ns Identificação de sistemas representados por equações a diferença 


j=4') e B(')=B(z'). Em seguida, determina-se x(t) pela 


Ale 
equação 


Alz Dt) -B(z Dult-d) 


% consequentemente, z(t) é obtida pela equação (5.37). 


b) Variável instrumental via medidas atrasadas 

Outro procedimento para seleção dos polinômios Ate'je Be!) 
da equação (5.37) é AE )=1 e B(2! )=z “» . Assim, o vetor z(t) pode ser 
reescrito como 


A0=[u(t-0) ... ult-/-na-nb) ; />d (529) 


que contém somente as entradas atrasadas (LJUNG; SÓDERSTRÓM, 
1983). A tabela (5.5) ilustra o código em Matlab do estimador da VI com 
medidas atrasadas. 


c) Variável instrumental refinada 

O método da variável instrumental refinada foi proposto por P. C. 
Young (1976), em que uma pré-filtragem nos dados do algoritmo é 
incluída. Portanto, o estimador de parâmetros é composto das seguintes 
equações básicas: 


t+) =) +Ke+D t+ D+ DÃO) 


P(t)z(t+1) 


NS ES TORTO PE 


(5.40) 


P(t+1)= á e as P(tz(t+ Dq (t+ o 


A + q (t+ DPCOalt+ 1) 

vd =T(z yo) 5 ql) =T(e Dto) 
Fr Te?) é o pré-filtro selecionado pelo usuário para evitar instabili- 
E ea e o vetor z(t) pode ser selecionado pelas equações 


Brtee 
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Tabela 5.5 - Programação em Mattab do estimador da Wi 


cê —— ” - os 
da variável instrumental 


» a uma planta de segunda ordem 


Quantas iterações ? VE 
% Entrada e perturbação 


1;p=1000*eye (4,4); % Condições iniciais 
;0;0;0]; zeta=[0;0;0;0): 


teta=[ 
for t= 
y(t)=0;x(t)=0;erro(t)=0; 
al(t)=teta(1);a2 (t)=teta(2);b0 (t)=teta (3) ;b] (t)=teta (4); 
end 
for t=6:nit % Iteração da simulação 
ylt)=1.5*y(t-1)-0.7*y (E-2) +u(t-1)+0. Stu (t=2) +. Er 
e(t)-e(t-1)+0.2*%e(t-2); 
fi=[-y(t-1);=y(t-2);u(t-1)sult=2)); 
zeta=[u(t-1)jult-2) ;u(t-3) u(t-4)]; 
erro(t)=y(t)-teta'*fi; 
k=p*zeta/ (lambda+fi'*p*zeta); 
teta=teta+k*erro(t); 
p=(p-k*£i'+p) /lambda; 


al(t)=teta(1);aZ(t)=teta(2);b0(t)=teta(3);bl (t)=teta(4): 
end 


telinit; 

subplot (221) ,plot (t,al (t)),ylabel('al'),x1 y a 
abel ('amostras 

subplot (222) plot (t,a2 (t)) ,ylabel ('a2!) , xlabel ('amostras")i 

subplot (223) plot (t,b0(t)) ylabel ('b0'), xlabel ('amostras')* 

subplot (224 lot (t (tj). ylabel ('b1'),xlabel ("amostras i 


Observação: 
“ Todas as seleções válidas da variável VI asseguram uma es” 


mativa não polarizada 
bilidade das estimatir assintótica de 6h. Entretanto, à 
ade 


pe vas é diferente cada escolha de 
Es que para uma dada quantidade de medidas, N 
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5.7.4 ALGORITMO DA VARIÁVEL INSTRUMENTAL 
SIMÉTRICA 


O algoritmo da variável instrumental apresentado pela equação 
(5.36) pode apresentar problemas em aplicações em tempo real causa- 
dos pelo mau condicionamento (definida negativa ou singular) da matriz 
de covariância P(t). Este problema pode ser contornado utilizando-se 
a fatoração UD de Bierman (1977), conforme discutido na seção (5.5). 
Para implementar a fatoração UD, é necessário introduzir uma nova 
forma para o algoritmo da variável instrumental, denominado algoritmo 
da variável instrumental simétrica (LJUNG; SÓDERSTRÓM, 1983) e repre- 
sentado pelas seguintes equações: 


due D= AO +KG+ yet+ =" (t+ 00) 


P(t)z(t+1) 


O e 
h+z (t+DP(tz(t+1) 


(5.41) 


A = sjro- Erg Dz" (t+ DP(t) 
À A +42" (t+ DP(DA(t +) 


3.8 ALGORITMO DE ESTIMAÇÃO DA MATRIZ 
ESTENDIDA 


] No projeto do estimador dos MQR, utilizou-se um modelo mate- 
mático para representar o processo do tipo CAR 


A(z Dy(D=2 “Blz Dult)+e(t) 


sm su a saída da planta está corrompida por um ruído do tipo 

PR anco (colorido), pode-se utilizar o algoritmo de estimação da 

ri estendida para evitar estimativas polarizadas (LJUNG, 1999). 
tindo um modelo média-móvel para o ruído e(t), isto é, 


D=C(z")w(t) (5.42) 
Burr 
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po so 
yaramétrico do processo caracteriza-se pel. - 
F pela Seguinte 


então O modelo 
função de transferência li 


near discreta: 
a(z y(D=2 “B(z Yu(t)+Clz Dw(t) 
y(0)=q"(1)0+w(t) (5.43 
43) 
onde os vetores de medidas e parâmetros são 
0(D=-y(t-1)=y(t=2)...— y(t-na) ult=d).. 
u(t-d-—nb) w(t-1)... w(t—nc)] 
BE(y=la, a; «asp by eis Ce C50] 


Assumindo que w(t) é mensurável, então a 
equai 9.43) é equi- 
valente à equação (5.5), ou seja, pads 


y(D=q"(D0+e(t) 


eos 1 

pia esa desconhecidos (polinômios estimados A(z?), B(z'), C(z')) 

pis Poa peeicrtog o problema do estimador dos mínimos qua- 
ep 4 ut) e w(t). Como na prática fw(t)) é um sinal não 

pois ser estimado (inferido) a partir do erro de predição, 


ni SO PO-F WE ps 
byiamente, i 
tis e considerar a parcela C(z:)w(t) como simples 


equação (5,43) util nc esereed dos polinômios (A(z"),B(z") da 
rs geradas são po) e o estimador dos MQR. Entretanto, as est 
ustra à Programação em M, s, à menos que C(z!) = 1. A tabela (5.8) 
dec mo ca fase] do estimador recursivo da matriz e 
ei cas estimador dos MOR (5.9) descreve o comportamento 
Presença de um ruído Mo 2 é de um parâmetro 
o. 
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ca 
Tabela 5.6 = Código em Matlab do estimador recursivo da matriz estendida 


“da matriz estendida 


put (! Quantas iterações 7 '); 
:nit 
nd>0.5 


if rar 
utij=l; 
else 
utij=-1; 
end 
end 
esu*0.1; 


1ambda=1;p=1000*eye (6,6) ;teta=[0;0;0;0;0;0); 


y(t)=0;w(t)=0;erro (t)=0; 
al(t)=teta(1);a2(t)=teta(2);b0(t)=teta (3);bl (t)=teta (4); 
cl(t)=teta(5);c2 (t)=teta (6); 


1.5*y(1-1)-0.7ty(t-2)+u(t-1)+0.5*u(t-2) te (t)-... 
e(t-1)+0.2*e(t-2); 
£i=[-y(t-1);-y(t-2);u(t-1)Gu(t-2);w(t=1);w(t=2)]; 
erro(t)=y (t)-teta'*fi; 
k=p*fi/(lambda+fi'*p*fi); 
teta-tetatkterro(t); 
p=(p-k*fi'*p) /lambda; 
wit)=y(t)-teta'*£i; 
al(t)=teta(1);a2Z(t)=teta(2);b0(t)=teta(3);bl (t)=teta (4); 
cl(t)=teta(5);c2(t)=teta(6); 

end 

t=l:nit; 

Subplot (221),plot (t,al(t)), ylabel ('al'),xlabel ('amostras"); 

Subplot (222) ,plot (t,a2(t)) , ylabel ('a2'),xlabel ('amostras"); 

Subplot (223), plot (t,b0(t))  ylabel ('b0'),xlabel ('amostras'); 


Subplot (224), plot (t,bl (t)) ylabel ('bl'),xlabel ('amostras!); 


rea TBlo 5.9 - Considere um sistema representado pela equação a dife 


y(t)=asy(t=1) + e(t) 
ênde a; = -0,6 e e(t) é uma sequência ruído colorido calculada pela se- 
Buinte expressão: 
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elo) = el) + 0.9e(t—1) 


do e(t) uma sequência ruído branco com média O (zero) e variância 
ra (5.15) ilustra à estimativa de a; pelo estimador dos MOR 
delo de primeira ordem da forma 


y(t)=á,y(t—1) 


sen 
1 (um). A figu 
admitindo-se o mo! 


Fiaura 5.15 - Comportamento polarizado na aplicação do estimador dos MOR 


estimativa do parâmetro a1 


observar o comportamento polarizado 


-0.798: é 
2), ou seja, o parâmetro estimado à: 


está incorreto, Portanto, o 
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5.9 PROBLEMAS 
1) Admitir o modelo matemático da forma 
v(t) = Bu(t) 


onde y(t) é a saída e u(t) é a entrada. Avaliar o parâmetro PB utilizando o 
estimador dos mínimos quadrados não recursivo para diferentes con- 
juntos de medidas, isto é, N = 3, 4, 5. As medidas de entrada e saída co- 
letadas são 


2) Identificar a função de transferência discreta de um processo gover- 
nado pela seguinte equação a diferenças: 


y(D+ay(t-1)=bou(t-1) 


utilizando o estimador dos mínimos quadrados (procedimento off-line). 
As medidas de entrada e saída são 


[a 
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| 0985 — 
-1.040 
1.256, 
14 | -0.212 0.966 


rar o processo discreto caracterizado pela equação 


7/8/9771 
EEE: 
-11/-08/-01| 0 


Estimar os parâmetros bo e b; pelo método dos mínimos quadra- 
dos não recursivo. 


3) Conside 


y(O= bou(+ byult— 1) 


4) Seja o processo discreto de primeira ordem CAR 
y(D=-a,y(t-1)+bou(t-1)+e(t) 


onde e(t) é um ruído branco com média O (zero). Uma experimentação, 


onde O <t <(N— 1), é realizada para a estima: do 
Os seguintes dados foram calculados: Sp or 


5 y(b=30 

Su(y=50 
Dyttei)y(d=1 
Ey(Ou(y=20 


o Dytt+ Du(t)= 36 
terminar pelo método 
mativas dos parâmetros oa ttlmas quadrados não recursivo as est 
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airutoS Identificação de sistemas repr 
capiruLo 5 ldentie 
E a ceê: iscreto de segunda ordem representado por 
tir o processo disc 
5) Admi 
0.52” Y(z) 
Gp(2z) a 


(1-0.52')(1-0.12") U(z) 


a) Obter a saída do processo v(t) para os seguintes sinais de entrada: 
* Sinal degrau: u(t) = 1; 
» Sinalrampa: u(t)=t; 
* Sinal senoidal: u(t) = sen(mt/5); 
b) Utilizando como modelo estimado a função de transferência discreta 
de primeira ordem 
boz 


G (= 
p(2) 1-a,z” 


determinar os parâmetros a; e bo pelo método dos mínimos quadrados 
(procedimento off-line) para as três entradas. Discutir a precisão dos re- 
sultados (admitir 20 medidas). 


8) Seja o sistema representado pelo seguinte modelo de segunda ordem: 


A(z"y(t=2"B(z")u(t)+e(t) 
onde 


A(z")=1-1.521 40.72? 
B(z1)=1+0.57* 
d=1 
Est “um sinal ruído branco com média O (zero) e desvio padrão 
“2eja u(t) uma sequência N(0,1) ou PRBS com amplitude 41, 


* Simular e identificar 
Md 9s parâmetros do processo com o estimador dos 
» Sie Para as seguintes condições: (na, nb) = (1, 0), (2, 1), (2,0), (2,1); 


Q tar os parâmetros estimados em função do nú 


N=1 00); mero de amostras 
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186 a = E 


) Analisar à variância e a função de covariância dos parâmetros esti. 
c 


mados pelo erro de estimação, isto é, 


e(b=y(t)- q"0 


7) Considerar um processo discreto de primeira ordem dado pela equa. 
ção a diferenças 
y(D=a,y(t-1)+bou(t—1)+e(t) 


onde os valores verdadeiros dos parâmetros são a; = 0.5, bo=1 eava- 

riância do ruído branco é 0.05. Avaliar a convergência dos parâmetros 

estimados sob as seguintes condições experimentais: 

a) u(t) = 0.0; 

b)u(t)=1.0; 

c) u(t) é um ruído branco com média O (zero) e variância 1 (um): 

d) u(t) = 0.1y(t); 

e) ut) = 0.1y(t) + r(t), onde r(t) é um ruído branco com média O (zero) é 
variância 1 (um), não correlacionado com o ruído do processo e(t); 

BD u(t)=0.1y(t-1); 

B) uít) = 0.1sinal(y(t)). 

O estimador dos mínimos quadrados recursivo deve ser inicializado com 

9s seguintes parâmetros: P(0) = 101», a(0) = 0, b0) = 1. 

8) Considerar o sistema de 


primeira orde; a inte 
equação a diferenças: m caracterizado pela seguin! 


Me y(D= Bleu 1) + lee) 


bo=2 para 0<t<920 
b=3 para 21<St<40; 
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Identificar os parâmetros do sistema via estimador dos MQR e comparar 

o desempenho dos estimadores com 

a) Fator de esquecimento: 2 = 1 e 0.9; 

b) Random walk: onde os elementos da matriz Q(t) são qu: = qu = qu = O 
equ=01el. 


9) Admitir um processo discreto AR dado por 
y(b=-a,yít—1) -a,ylt —2) + e(t) 


onde os parâmetros são a; = —1.5 e a; = 0.56. O ruído branco e(t) apresenta 
média O (zero) e variância 1 (um). Admitir a;(0) = a(0) = O e P(O) = 10 lx 
para 300 iterações. Comparar o desempenho do estimador dos MQR-pa- 
drão com o algoritmo da aproximação estocástica de modo que 


a) (a, /8) = (0.1, 1) e (0.01, 1) (sequência ganho com aprendizado sequen- 
cial); 

b) y=0.005 e 0.01 (sequência ganho escalar); 

c) yajustado pelo vetor total de medidas. 


10) Considerar o sistema discreto de segunda ordem CARMA caracteri- 
zado pela seguinte equação a diferenças: 


y(D-1.5y(t-1)+0.7y(t-2)=u(t—1)+0.5u(t—2) + 
e(t)-e(t-1)+0.2e(t— 2) 


o u(t) é um sinal do tipo PRBS e S/R = 10. Comparar o desempenho 
estimadores da variável instrumental e da matriz estendida, onde 


3) Para o algoritmo da variável instrumental o modelo a ser estimado é 
VD +á,y(t=1) +á,y(t=2)=b,u(t=1)+b,u(t—2) 
Avariável instrumental é selecionada por 


=) =[=x(t-1) =x(t-2)u(t- 1) u(t— 2)" 
x(t)=2"()0(t-1) 
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ao) = 0, P(O) = 1000Lsxay A = 0.99 e N = 100, 500, 2000 iteraçã 

tilizar O(t — 1) em lugar de O(t) é uma tentativa de a 
dade e convergência do algoritmo; a 
dida o modelo estimado tem a forma 


onde 
A razão para U 
lhorar a estabili 
b) Para o algoritmo da matriz esteni 


y()+á,y(t-1)+à,y(t=2) =bu(t-1)+b,u(t-2)+ 
dt) + e(t=1)+,e(t—2) 


onde 0(0) = 0, P(O) = 1000Iaxé, À = 0.99 é N = 100, 500, 2000 iterações 
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CONSIDERAÇÕES 
FINAIS SOBRE 
IDENTIFICAÇÃO 


6.1 FUNDAMENTOS PARA IDENTIFICAÇÃO 
VIA RELÉ 


Na atualidade, em função de uma grande parte dos controladores 
manufaturados de processos industriais apresentar um sistema de 
sintonia denominado auto-tuning, é inevitável uma breve descrição da 
metodologia do relé no contexto da estimação de modelos matemáticos 
de ordem reduzida a partir da resposta em frequência (HANG et al. 1993; 

ÓM; WITTENMARK, 1995; ALMEIDA; COELHO, 2002; VISIOLI, 
2006; YU, 2006). 

* Os experimentos com relé na malha de realimentação, com o pro- 

pg de identificação de processos, tornaram-se populares a partir do 

| pbalho de Astróm e Hágglund (1984). Este método foi utilizado pa- 
| determinar o ganho crítico e a frequência crítica e, por consequência, 


] jd 
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automatizar O método 


i icho! 
to por Ziegler e Nic 
da ed linearidade atra 


de oscilação de projeto do controlador PIp ” 
Is (1942). A abordagem baseia-se na moda E 
vés de sua função descritiva e na interpres q 
ermos do diagrama de Nyquist para a obtenção de informa 
cia do processo. À identificação do processo é Es 
ação em frequência da função de transferência y 

o 


perimento de malha fechada. A figura (6.1) ilustra 


em t 
em frequên 
partir da estim 


processo num ex) 
elemento não linear (relé) realimentando um sistema de controle para 
0 


propósito da identificação paramétrica (relay-feedback experiment) 


Figura 6.1 - Realimentação do processo através do relé 


ge fRsD= = 


; = ão do processo e da especificação da não linearidade, deter- 
po à Re relevantes (amplitude e frequência de osc- 
o so à estimação da função de transferência do pro- 
tipos de relés: a) o relé ja). a à estimação de G(ja), utilizam-se dois 
Popeda à e e histerese é utilizado para estimar a função 

quência de cruzamento, ou seja, na parte nega 


tiva do eixo real ' 

para estimar a dedo de Gio); b) O relé com histerese é empregado 
seleção do relé é determinada Pç cia em diferentes frequências. à | 
Planta, da precisão do modelo mai ção das condições operacionais | 


se [e(t) > 0 


se fe(t) « o RtÃo “u(t) = 4d 


É então u(t) = -d 
o rel 
do com histerese, figura (6.2b), pode ser mo 


vendo o co tempo i 
Portamento da beber pone regras linguísticas desce 
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sobre identificação 


se [le(t)| > E 8 e(t) > 0] então u(t) = +d 


se [le(t)| > & & e(b) < 0) então u(t) = -d 
se [le(D)| < e &eu(t= 1) = +d) então u(t) = +d 
se [le(t)| < E &u(t - 1) = -d] então u(t) = —d 


nde de e são definidos conforme a figura (6.2b). 
o 


Dee 


| Figura 6.2 - [a] Relé sem his se; [b) Relé com his 


uít) u(t) 
d 


ett) 


(a) (b) 


Figura 6.3 - Saída do relé u/t/, sinal quadrado; 
saída do processo yit) característica senoidal 


e rejé sem histerese relé com histerese 
My 15 


Ts 
5 
Bs 
g 
8 
8 
s 
B 
g 
8 
g 
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DEM DINÂMICOS 


O LINEAnge 


Ea? ye u(t) são mostrados na figura (6.3). A saída do 
Os sinais vt e sponde a uma onda quadrada, Com 
Com os 


É -ontrole, corre 
ult), variável de cont 3a) e (6.3b), como entrada para q p 


s. figuras (6 
ad jo relé, figur 5 
sinais de na se que a resposta do processo de malha fech 
sidero x 


cesso, e cons tes de baixa frequência, a saída osci 
2 : las componen es es oscil, 
seja eia aÃ (6.3a) e (6.3b). 

forma seno dal, 


relé 


ro- 
ada 
a de 


Observações: 

Na prática, a histerese é implementada no relé, conforme à 
figura (6.1), uma vez que, na presença de ruído e sem a his. 
terese na malha, o relé pode chavear nos estados liga/desliga, 
isto é, apresentar um comportamento oscilatório randômico; 


* Em aplicações práticas, a largura da histerese é selecionada 
com base no nível do ruído, por conveniência, duas vezes maior 
do que a amplitude do ruído (HANG et al., 1993); 


* O método clássico da resposta em frequência de Ziegler/ 
Nichols tem algumas vantagens. A técnica baseia-se em um 
simples experimento e procura identificar a frequência de os- 
dilação (ultimate frequency) da saída do processo. Entretanto, é 
difícil automatizar ou executar este experimento de tal for 
Vi sadia ie ses oscilação mantenha-se sob controle 
E " » O funcionamento do processo próximo € 
cnh tg fazendo com que seja ni ár 

: uxos normais de operação de uma planta &” 


WITTENMARI 20 é OM; 
MARK 2001) credível (OGATA, 2010; ÁSTR 


Vários r 
Pontos sobre a curva de Nyquist podem ser determi 


n 
ções entre Fr à da experimentação com diferentes 

Siclo-limite em ai né Nalmente, controla-se a amplitude 
amplitude do relé, vel desejado pela seleção aproprial 
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6.2 ESTIMAÇÃO DE MODELOS PONTO A PONTO 


A função descritiva ou função descritiva senoidal de um elemento 
ão linear é definida como a relação complexa entre a componente har- 
o ? 4 
o anica fundamental do sinal de saída e do sinal de entrada, isto é, 
mo 
x 
N=-— £FASE 
U 


onde N é a função descritiva, U é a amplitude da componente fundamen- 
tal do sinal de entrada e Y é a amplitude da componente fundamental do 
sinal de saída (OGATA, 2010). 
Considerando o relé sem histerese e com histerese, têm-se as se- 
guintes equações relatando as funções descritivas: 
1 ra 1 T ZA rE 


* ut aa O Uta 


“NG) 4d 


A partir da modelagem do relé por função descritiva e da operação do 
sistema sob o controle do relé, pode-se determinar a função de transfe- 
rência do processo por 
! a ma 
SG, (j0)=-— > G,G0)=-— 
plo) N(a) p(jo) 4d 


ja ME 
acne 


onde a é a amplitude de oscilação do sinal na saída do processo e «é a 
frequência de oscilação medida. 

Pela equação (6.1), pode-se estimar a função de transferência do 
Processo na frequência de cruzamento utilizando-se o relé sem histe- 
ge À equação (6.1) permite também estimar a função de transferência 

Processo em diferentes frequências utilizando-se um relé com his- 
pe € ajustando-se diferentes valores para o parâmetro 6; As inter- 
“t0es, no plano. G;(s), dos lugares geométricos do recíproco inverso da 
função descritiva do relé com o lugar geométrico de G;(j0)) estabelecem 


(6.1) 


Bone 
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Vi Tal 
o correspondentes às frequências da parte a 
tuadas no terceiro quadrante do plano Gs) ne 
das as frequências que, em geral, são de A este 
sistemas de controle e que ico 
és de experimentos com o relé, a 
a 


pontos de operaçã 
grama de Nyquist si 
quadrante, estão situa 
se para análise é projeto de 
completamente definidas atravi 


figura (6.4). 
Figura 6.4 - Interseção dos lugares geométricos do recíproco inverso da função de 
do relé com histerese la] e sem histerese [b] com o lugar geométrico de EM g 
» IMG, (ju) 4 
AÍ É ImS,(jw) 


ReG (jo) 


ç (b) 


Observações: 


as representação do si 

drada) pelo pila de saída do relé (que é uma onda qua- 
descritiva) é uma aproxi da série de Fourier (ou pela função 
temas de ordem Es ção que, para muitos processos (sis: 
promete o de: head atrasos significativos, etc.), com 
1997a, 1997). O controlador sintonizado (WANG eta. 
dm ponto da fu peso possibilita a identificação de apa 
pode ser insuficiente para pa em frequência do sistema, o que 
malha fe dO em falha na Screver satisfatoriamente o proce. 

fechada indesejável; sintonia do controlador e dinâmica de 


elevado ní 
vel de ruído Fe caça utilizada em sistemas co” 
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jáveis, uma vez que o ponto de detecção do chaveamento do 
relé é corrompido pelo ruído presente no sinal; 

Para a estimação do ponto de interseção da função de trans- 
ferência com o eixo imaginário negativo, utiliza-se um inte- 
grador com o relé na malha de realimentação (sem histerese). 
Neste caso, a função de transferência é calculada por 


Eq) == [520 
= (6.2) 


A tabela (6.1) ilustra a programação em Matlab da técnica do relé 
com histerese na modelagem ponto a ponto do processo de segunda 
ordem do exemplo (6.1). 


Exemplo 6.1 — Seja a seguinte planta contínua de segunda ordem com 
atraso de transporte: 


G,(s) 


e -28 


“(255+1)6.755+1) 


Identificar o processo em diferentes frequências utilizando a téc- 
nica do relé onde d = 0.5 e £= 0.0307, 0.0474, 0.0870, 0.1632. 


Note que os erros estimados da frequência e da parte real man- 


a 
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ica do relé 
Tabeta 6.1 - Código em Mattabdatécnicadorelé 
a artes real E imaginária 


% Planta de 
clear all;close allicic; 
Tamostra= nptos=100:d=0.57 
for t=1:7, 
ut) e-dye(t)=05y (t)=0: tempo (t)=E 
end; 
for t=5Sinptos, 
y(t)=1.436%y(t-1)-0.5134*y (6-2) +. ++ 
0.04286*u (t-2)+0.03431%u (t-3) 5 
el(t)==yít): 
if ((abs(e(t))>eps) & (e(t)>0)) u(t)=+d; end; 
if ((abs(e(t))>eps) & (e(t)<0)) u(t)=-d; end; 
1£ ((abs(e(t))<eps) & (u(t-1)==+d)) u(t)=+d; end; 
if ((abs(e(t))<eps) & (u(t-1)==-d)) u(t)=-d; end; 
if (e(t)==+eps) u(t+1)=td; end; 
if (e(t)==-eps) u(t+1)=-d; end; 
tempo (t)=t*Tamostra; 
end; 
ch=0;i=1; 
for k=4:nptos, 
sl=u(k);s2=u(k-1); 
if sl-=s2 
ch=ch+1; 
i=i+l; 
end 
end 
per=( (nptos+ “9% 
ns 4)/ch)*2 Tamostra;omegas (2*pi) / (per); 
for t=l:nptos, % Valor de pico 


i£ y(t)>=arm arm= 
end; 


a=arm; 
realG=- (pitsgrt (a"2-epas 
imagG=- (piteps)/ 4R) a (44d); 
saida=[u;y];plot (tempo, saida) 
; 


eps=0.0870; 


*Tamostra; 


% Experimentação com o relé 


e 


% Período 


Y(t); end; 


% Partes real e imaginária 
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6.2.1 ESTIMAÇÃO NÃO RECURSIVA PELA 
RESPOSTA EM FREQUÊNCIA 


O conhecimento da função de transferência de uma planta é 
necessário em diversas aplicações (predição, controle, supervisão, diag- 
nóstico, etc). A seguir, apresenta-se um método não iterativo para 
estimação em frequência de um modelo contínuo estável de segunda 
ordem com atraso de transporte com base nas partes real e imaginária 
do sistema, isto é, magnitude e fase (WANG et al., 1997b). 

Considere o seguinte modelo matemático da planta: 


Ls 
e 


G()=—— 
E as? +bs+c 

que pode representar processos monotônicos ou oscilatórios. As equa- 

ções de magnitude e fase são dadas por 


1 


CU === 
E, jo) Vc-am?)? + (bo)? 


jl==io= Eae 
arglG, Go)]= Lo arca DO) 


Se a resposta em frequência do processo G;(jcu), i = 1, 2, ..., M está 
disponível a priori, por exemplo, calculada pela técnica da estimação de 
modelos ponto a ponto, conforme apresentado na seção (6.2), então é 
possível reescrever a equação da magnitude na forma clássica do esti- 
mador dos mínimos quadrados não recursivo, ou seja, 


ade P=40 
Pesto, of ai 1] 0 ab 
' 1 
5 /8,G0)! |; q= a | ; 05] 0, |=| (b?-200) 
: E) : 0 c 
De] 
E, Goo Bones 
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s do modelo a, b, 6, L, isto é 
e, assim, as parâme nodelo a, D, €, L, 1s 5 
im, i ificar 'àr tros d 

, assim, identificar os f 


se J0,0, |; cam?) ) 


711 0 nf bo 3 
al 1 E 
N 0, = 0,0, |; Lo= arglG, (0)] Eira ( 
4) 


| 


Observação: 

* Aplicando-se o método de estimação ao processo de segund; 
ordem do exemplo (6.1), onde d = 0.5 e & = 0.1931, 0.1632, 
0.1355, 0.1101, 0.0870, obtêm-se os seguintes valores para- 
métricos: a = 11.2562, b = 7.7003, c = 1.3542, L = 2.3088. Logo, 
a técnica mostra-se atrativa para a modelagem de sistemas 
através do experimento via relé e pode também ser aplicada 
na sintonia automática de controladores de processos. 


6.3 ESTIMAÇÃO DE MODELOS DE PRIMEIRA 
E DE SEGUNDA ORDEM 


Recentemente, o método de identificação através do relé vem 


sendo pingndo Fa a Ei de modelos matemáticos de baixa 
ordem à atraso de transporte i e 
industriais (WANG et al,, 19973). Um Para aplicações em process 


Control), entre outros (WR, é Smith, GPC (Generalized Predictive 

ano Pp VISTOLA 206? ZARROP, 1991; HANG et al, 1993 
a prática, a ] 

por 1 cia dois dos Processos industriais pode ser 

dem com atraso de tanepor anda Primeira ou de segunda or 

19974; ALMEIDA, 2002), Cop AO RARI; ZAPIROU, 1989; WANG et al: 

está representado por "ando que o modelo estimado do pro” 
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Ke 


sas ts +1 


(6.3) 


corelé é implementado conforme a figura (6.5a), então expressões para 


o at 


(1997a). 


Figura 6.5 - (a) Relé não simétrico 


[b] Saída do relé londa quadradal e saída do proc: 


u(t) 4 


do + d 


et) 


o 5 10 15 20 25 30 35 40 
tempo (5) 


raso de transporte (9), ganho estático do processo (K,) e constante 
de tempo (1) podem ser obtidas, conforme apresentado em Wang et al. 
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Seja o ganho estático do processo calculado por 
Lat Ts 
EredE 
K k TO =6,(0) 4 
ll “u(t)dt 
o 


onde Tu e Tu estão definidos na figura (6.5). Assim, a constante de 
tempo pode ser obtida por 
Tu 
a E 
i ta E +d, JK, e” —(d, —do)K, 80 +E 
(dy +d,)K, (E, +8) 


(6.5) 


onde ds, dh, d», &, & são definidos na figura (6.5a) e o atraso de trans- 
porte resultante é 


' mp So | 
1 


A td, -doK, x 


— A figura (66) e a tabela (6.2) ilustram a programação em Matlab/ 
Simulink da técnica do relé com histerese na estimação de modelos de 
primeira ordem aplicados a um processo de segunda ordem. 


Figura 6.6 - Diagrama em Simulink da modelagem via relé 


o] E] 
dE 


-a 


pd. 


sinal de 
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aan e A 
Tabela 6.2 = Código em Matlab: relé não simétrico 

% Método do relé para identific de um modelo 
» de transporte, isto é, 


de primeira ordem com atr 
iotexp(-tetal*s)/(taults+1) 
ur e yr são obtidos através do Simulink 
nd) 


4 


vetor 
di=0.5;d2=0.5/eps0=2;epss= 
O; Tamostra=0.1; 
nptosl=nptos/Tamostra; kont=0; 
for t=2:nptosi, 
if (ur(t)-=ur(t-1)) 
kont=kont+l; ch (kont)=t; 
end; 
end; 


Tul=(ch(6)-ch(5)) *Tamostra; Tu2= (ch(5)-ch(4)) *Tamostra; 
auxi=ch(3) ;aux2=ch(5) ;i=0; 
for t=auxl:aux2; 

l=i+; 


E(t);ui(i)=ur(t); 
*Tamostra; 


SH (10 yil+tyi 0]).*([ti 0J-[0 til); 
um(al(1,2:1ength(yi))); 
a2=0.5*([0 uíi]+[ui 0]).*([ti 0]-[0 tãj); 
aZ=sum(a2 (1,2:length (ui))); 
K=al/a2 % Ganho 
arm=do; 
for t=aux] :aux2, 
pis Yr(t)>=arm arm=yr(t); end; 
Ausarm; 
arm=do; 
for t-auxl :aux2, 
Pri Yr(t)<=arm arm=yr (t); end; 
Adsarm, 
tetasjog ” * Pico negativo e Atraso 
*Le(dLrdz) datado) *K-epsteps0) /((dl=do) *Ko+Ad)); 
“2 (dosao Krfexp (teta) - (dl-dO) *Kpteps-epsO; 

2 V*k-eps-eps0; 

“Tul/1og x1/x2) 


( 
Eetaleta * Cte, de tempo 
siga % Atraso de transport 


% Pico positivo 
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Outro método disponível na liter 
para estimação do modelo matemático 
baseia-se nas seguintes equações 


Observações: 


ratura de controle de processos, 
FOPDT, conforme equação (6, 3) 
(HANG et al. 1993): 


ca 2 
Va JK, 1 (6.7) 


ue 2n1, 
=-2| q-arctan——— (6.8) 
0, el? ar T ) 


u 


O ganho estático K, pode ser estimado a partir das medidas de 
entrada e saída em regime permanente para uma mudança 
degrau na entrada da planta; 
Se o relé sem histerese é empregado, o ganho final K, é calcu- 
lado por 

o ses 


ra 


enquanto que, se o relé é implementado com histerese, então 
utiliza-se 


pese dds. 
Coma-s' 


Os parâmetros K, e T, = Tur Ta (período final) são obtidos do 


experimento com q relé, 
equações (6.7) e (6.8), SE 8 À adam ser calculados PM 


à O, aa 
pregar uma Era à calibração do controlador PID pode em” 
Ta = T/8 (original) es fórmulas: K = 0.6, Ti = 1/2 

| K = 0.33K Ti = T/2, Ty = T/3 (algum 
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sobressinal); K = 0.2K Ti = Tu, Tá = To/3 (sem sobressinal) 
(YU, 2006); 

» Para o processo modelado pela função de transferência de pri- 
meira ordem, equação (6.3), as regras empíricas de sintonia 
de Ziegler/Nichols ou Cohen/Coon podem ser aplicadas na pa- 
rametrização do controle PID. As tabelas (6.3) e (6.4) ilustram 
a sintonia (tuning) do controlador PID (BROSILOW; JOSEPH, 
2002; ALMEIDA, 2002; VISIOLI, 2006; MARIO, 2008). 


Tabela 6.3 - Parâmetros PID com base na sintonia de Zregler/Nichols 
Ape é epic sed 1<0/4 «nl 


(UK D(/0) 
(0.9/K (11/03) 
PID | (.2/K)1)/09) 


30+3(0,/7,) 
*9+20(0,/7,) 


32+6(0,/t,) o, 


*13+8(0,/7,) ep tr) 


atraso de 


Considerando agora um modelo estimado de segunda ordem com 


para a planta, então tem-se 
perto 


(6.9) 
(tos + 1)? 


G,(9)= 
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ransporte /) estão relacionados 


so atraso de t es 
esa lo final T, pelas Seguin. 


de tempo 7 : 
À constante F ENE final Ky e períod 


com o ganho estático Kos ka 
tes equações (HANG; CAO, 1996): 


A Ee 
=> KK,-1 (6.10) 
t Fi) ' 


g,=l x-2artandE | Tefa-2arctan( ER, —D) (611) 
“ 2a T, 2x 


u 


onde K, e T. são obtidos do experimento com o relé e 7 e 6 podem ser 
calculados a partir das equações (6.10) e (6.11), respectivamente (HANG 
etal., 1993). 


Observações: 

* A aplicação dos modelos de planta estimados com ordem 
reduzida (primeira e segunda ordens) não se limita ao projeto 
de controladores PID, podendo também ser utilizados na sín- 
tese de controladores clássicos ou avançados; 


* Os modelos das equações (6.3) e (6.9), e as correspondentes 


fórmulas de sintonias, são comuns em processos industriais 
adaptativos (auto-tuning, selftuning); 


de controle está nom 


está no estado 4 
determinados E poripio constante, os parâmetros PID são 


rmazenados E 8 
dos automaticamente (os valore: 
Parâmetros PID pr atualizados com o sinal do exper 
on-demand). ntrole também é denominada tuning 

je Po do processo mudam, então O 
2003; MARIO, 2008). P9de Ser ressintonizado (VANDOREN, 
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E at, o] PD [uno [o 
= y e +) proseso | att, 


6.4 APLICAÇÕES VIA EXPERIMENTAÇÃO 
E SIMULAÇÃO 


Exemplo 6.2 — Um exemplo da aplicação prática do método do relé em 
um processo denominado fan-and-plate, desenvolvido no Laboratório de 
Controle e Automação do Departamento de Automação e Sistemas da 
5 Universidade Federal de Santa Catarina (LCA/DAS/UFSC), está ilustrado 
na figura (6.8). Informações adicionais estão disponíveis em Coelho et al., 
2001). A seguir, faz-se uma breve descrição do processo (protótipo físi- 
co) e aplica-se o método do relé na etapa de identificação. 

Um motor DC acelera o jato de ar através de um túnel de caracte- 
rística afunilada. Na extremidade do túnel, encontra-se posicionada 
uma placa retangular. O objetivo de controle de malha fechada é regular 
à deflexão angular da placa. O tamanho do túnel e as entradas laterais 
de ar podem ser regulados, provocando perturbações de carga. Pertur- 
Fi elétricas também podem ser adicionadas ao sistema. O processo 
Deda nedeada por fase não mínima, atraso de transporte e fluxo tur- 

to e ressonante. Assim, pode-se validar técnicas de identificação 


“ob difíceis condições de operação em um sistema complexo. 


= 


ss Es 
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co do proce 
g - Diagrama esquemát 


E 
instrumentação 


late, o resultado do experimento com o 


Para o processo fan-and-pl Ret 
relé e o diagrama de Nyquist referente ao modelo de primeira ordem 


com atraso de transporte são apresentados na figura (6.9). Neste experi- 
mento, utiliza-se um relé com histerese segundo as seguintes especifi- 
cações: d; = d» = 1, & = 0 e &= 0.3. Um período de amostragem de 0.2 s 
é utilizado na aquisição das medidas e o processo é fixado no meio da 
faixa de operação (do = 2.5 volts). 

Como resultado da aplicação das equações (6.7) e (6.8), o modelo 
matemático estimado de primeira ordem com atraso de transporte para 
o processo fan-and-plate é expresso por 


a -0,21s 
89= 0.822e 
0.12s+1 


E ' 
para o Eine ese estimado de baixa ordem tenha sido obtido 
xo nd ds “Plate, ele não representa adequadamente o pro 
9s níveis de funcionamento. O modelo matemático de 


primeira ordem representa xima: 
modelo do processo na faixa Pá ei E rija ese 


cedimento de identificação. Snamento do processo durante 0 pro” 
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Figura 6.9 = Aplicação do relé no pro pl 
$ | 
| Ns ída di | 
05h dotes fanandplate -| 
E —80 E) E 20 TOS 
tempo (s) 
(a) 
44 ImG,(s) 
> 
-0.8 ReG,(s) 
W=8.62 radis 
U=4.48 radis 
w 
(b) 
iação do método do 
Exemplo 6.3 - Como exemplo de simulação asi gre mi de baixa or- 
telé (não simétrico com histerese), através ep ado É 
+ Considere um processo contínuo des 
-28 
aqua 
Go(8) (Bs +17 
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e os parâmetros d 


d;= 0.5. A partir das equaçõe 


primeira ordem para O processo € 


idas 


ImG,s 


o relé ajustados com 
es (6.4), (6.5) é (6. 


6,(9)= 


no: & = 


273269 


3.4094s 41 


do=2,d=054 
6), o modelo estim; q 
“timado de 


40 
tempo (s) 


50 


so 


70 
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A figura (6.10) apresenta º processo sob o controle do relé, 
enquanto que a figura (6 10b) apresenta as curvas de Nyquist real e esti- 
mada para o processo. Observa se pela figura (6.10b) que, do ponto de 
vista do diagrama de Nyquist, resultados satisfatórios são obtidos nas 
baixas, médias e altas frequências, resultando, portanto, uma estimação 


adequada. 


Figura 6.11 - Simulação de malha fechada do processo de fase não minim 
je fase não mínima 


s 
v 
5 J 
õ 
E) 1 
o 
s J 
E 
J 
so 
v 
£ 
5 
8 
0 10 20 30 40 50 [5] 


tempo (s) 


Exemplo 6.4 — Os sistemas de fase não mínima representam um desafio 
à engenheiro de controle de processos devido à resposta inversa inicial. 
E se citar como exemplo de sistema de fase não mínima o controle 

a de água em caldeira geradora de vapor. A seguir, realiza-se um 
“Perimento onde se observa a dinâmica temporal imposta pelo relé e, 


em Seguida, a resposta para mudanças de referência do sistema de 


e ole PID ressintonizado. A figura (6.11) mostra o resultado temporal 
a A numérica PID auto-tuning implementada com o processo 
“se não mínima modelado por as 
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A tabela (6.5) ilustra o código de programação em Matlab na 
sintonia e controle realimentado que obedece à seguinte parametri- 
zação: Ku = 1.4895, Tu = 1.1 51 (dm Amin) = (4, 0), Ti = 0.055, y,=2, 
sintonia PID inicial (K, To To) = (0.2, 0.75, 1) e sintonia PID final (K, T, 


Tabela 6.5 - Código de programação em Matlab do PID auto-turing 
E Modelagem da função de Fransterência: Relé sem histerese 
3 Controlader PID Ideal Incremental 
% Função de transferência continua 
5 (-0.228 + 1) 
g E 
t (0.225 + 1)º2 


$ 

clear alliclose all;clc; 

G(1:4) =0; du (1:4)=0; umin=0; umax=100;Ts=0 .05;nptos=1100; 

Kc=0.2;Ti=0.75;Td=l; 

qgo=Kc* (1+ (Ts/Ti)+ (Td/Ts) ) pal=-Ke+ (1+2* (Td/Ts)) ;q2=Kc* (Td/Ts); 

ref=2;dnax=ref+ref; dminsref-ref;d=dmax-dmin) 5 

y(1:4)=0;e(1:4)=0; 

for t=1:4 tempo(t)=t*Ts; end; 

4 Discretização da planta 

nps=[-0,22 1) ;dps=conv([0.22 1], [0.22 11]); 

[npz, dpz!=c2dm (nps dps, Ts,'zoh') ; 

yr (1:700)=ref;yr (101:900)=2*ref; yr (901:1100)=ref; 

4 

for tsfnptos 

£ 

if (1<=400) | (t>=601) 
ytt)=-dpz (2) *y(t=1) = av (t- = 

E )=dpz (3) *y(t-2) +npa (2) tu tt=1) te 

elt)=yr(t)-y(t)s 
ulc)=u(t-1) +q0te (t) tglte(t-1) +q2te (t-2) 7 


elseif (t2=401) 4 (t<=600) 4% Planta com relé 


(r)=-dpz (2) *y(t-1)-: 
piRiean AMRS dE CER (3) dy Ei feio CREA à 


e(t) = yrtt) - yit); 
if el >=0) ulti=dmax; end; 


E Ot) Cori cd ai ea a o da 
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JEle(t-l)< 0) ult)=dmin; end; 7 


tempo (t)=t*Ts; 


endt 

(1 1£ t==600 
1] kont=0; % Calcular período 
» for t=401:600, 
pé if(u(t)-=u(t-1)) 
À kont=kont+1;ch (kont)=t; 

end; 

end; 


Tul=(ch(7)-ch(6))*Ts;Tu2=(ch(8)-ch(7))*Ts; 
Tu=Tul+Tu2; 


Wu=(2*pi)/Tu 
auxl=ch (5) ;aux2=ch (7) ; 
armeref; % Calcular valor de pico positivo 


for t=auxl:aux2, 
if y(t)>=arm arme=y(t); end; 


end; 
Au = arm; 
arm = ref; % Calcular valor de pico negativo 


for t=auxl:aux2, 
if vy(t)<arm arm=y(t); end; 
end; 
Ad=arm; 
a=abs (Au) -abs (Ad) ; 
Ku=(4*d) / (pita) 


Kp=1; % Calcular parâmetros do modelo 
Kc=0.2*Ku % PID Ziegler-Nichols (sem sobressinal) 
Ti=0.5*Tu 

Td=0,125*Tu 


q0=Kc* (1+(Ts/Ti)+(Td/Ts)) ;ql=-Kc* (1+2* (Td/Ts)) 5 
q2=Kc* (Td/Ts) ; 
end 
$ 
end 
t 


tempos0:Ts 


(nptos-1)*Ts; 


Sacos [u; 
Subplor (2,1,1) plot (tempo, yt, !== 
| k' tempo, y, et! Pi inemident: 2) , Vlabel ('satda e referência"), 
*label ("tempo (8) '), 
grid on; 
Subplot (2,1,2),plot (tempo,u, 'b', 'Linewidth',2), 


Ylabel ('controle"), xlabel (* (s) ')y grid oni. E 
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6.5 PROBLEMAS 


içô ão e os sistem; 
1) Considerar as seguintes condições de operação as de con- 


trole não lineares: 


a) Relé sem histerese com amplitude d=1ea planta de terceira ordem 


caracterizada por 
SE » 4 
rs) s(s+1)(s+2) 


b) Relé com histerese de parâmetros d=5,E=1 ea planta de primeira 
ordem representada por 


(=) 2541 


Para cada sistema, determinar a amplitude e a frequência do ciclo-limi- 


te da saída. Verificar os resultados com o ambiente de simulação Matlab/ 
Simulink. 


2) Obter as constantes K, e Tu a partir da figura (6.3). Comparar e co- 
mentar os valores obtidos. 


3) Utilizar os valores de K, e Tu encontrados da fi; - 

d E igura (6.3) e, const 
derando o ganho estático do processo K, = 1, determinar o modelo de 
primeira ordem com atraso de transporte da equação (6.3). 


4) Seja o processo de segunda ordem descrito pela seguinte função: 


nte 205 
há “(s+1) 


Se o relé sem histerese tem d = 0, 
se obtêm K,= 4.2, Tu = 3.3€ 0 modao pe eoS: Mostrar por simulação que 


modelo FOPDT tem a forma 
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o Um processo de segunda e der piaspomdançã é comum em pro- 
: jais. O seguinte modelo é simulado; 
industriais. 
cessos 
il 
G(s) 


PE (6+1)(5.175+1) 


m ganhos iniciais do PID de valores K. = 1.2, T; = 0.75, Ta = 0.75. 
o que a resposta temporal de malha fechada obedece à figura 
(6.12) para à parametrização do algoritmo do controlador PID ideal 

1, ou seja, dna = 6, dmin = 2, Jr = 4, Ti = 0.55, K = 7.6912, 


autoajustável 
T,=6s, com a sintonia final obtida para K. = 1.5382, T;= 3, Ta= 0.75. 


Figura 6.12 - Simulação de malha fechada do processo de segunda ordem 


s 
é 
z 
v 
g 
E) 
| 
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